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SECTION  I. 


INTRODUCTION 


The  need  for  analyzing  images  of  the  earth  taken  from 
high  and  low  altitude  aircraft  is  greatly  increasing.  The 
purposes  of  the  photographs  are  diverse  and  include: 

.  Navigational  aid 
.  Geographic  map  making 
.  Natural  resources  analysis 
.  Weather  prediction 

Typically,  a  sequence  of  images  is  taken  at  different  times 
with  different  sensors.  A  basic  problem  in  the  analysis  of  these 
images  arises  when  it  is  desired  to  locate  the  same  region  in 
a  pair  of  dissimilar  images  of  the  same  field  of  view.  For 
example,  one  may  wish  to  locate  a  feature  in  both  optical  and 
radar  images  of  a  scene  in  order  to  verify  the  accuracy  of  the 
radar  sensor  or  locate  a  feature  in  two  optical  images  taken 
under  different  weather  conditions  in  order  to  study  the  changes. 
This  general  problem  is  called  map  matching.  Special  cases  of 
this  problem  are  image  registration,  change  detection,  stereometric 
analysis,  and  multispectral  image  analysis. 

Several  methods  including  photographic  techniques  such  as 
image  subtraction,  optical  techniques  such  as  optical  correlation 
and  digital  computer  techniques  may  be  used  to  approach  the 
map  matching  problem.  In  this  report  only  computer  imaging 
techniques  will  be  considered.  The  versatility,  reliability,  and 
economical  feasibility  of  these  techniques  have  produced  the  trend 
of  these  methods  for  all  application  areas. 


Several  noteworthy  approaches  to  the  map  matching  pro¬ 
blem  have  been  proposed.  These  include  the  standard  correlation 
detector  as  well  as  highly  efficient  sequential  techniques.  These 
methods  and  techniques  together  with  an  analysis  of  their  perfor¬ 
mance  are  described  in  the  next  section.  The  results  of  this 
evaluation  indicate  that  although  some  of  the  methods  are  partially 
successful,  further  research  is  needed  to  accomplish  the  difficult 
task  of  map  matching  in  an  efficient  manner. 

The  type  imagery  encountered  in  map  matching  is  varied  and 
may  include  optical,  radar,  radiometric  or  multispectral  images. 

A  typical  optical  image  is  shown  in  Figure  1,  and  the  correspond¬ 
ing  side  looking  radar  image  is  shown  in  Figure  2.  Note  that 
although  each  of  these  images  shows  significant  scene  detail,  the 
two  sensors  respond  differently  to  the  same  feature.  Most  noti¬ 
ceable  is  the  intensity  reversal  of  roads  which  appear  white  in 
the  optical  and  dark  in  the  radar  image.  Upon  careful  study  one 
may  also  note  a  nonlinear  geometric  distortion  between  the  two 
images  as  well  as  scale  changes.  The  appearance  of  bright  "glint" 
features  in  the  radar  image  may  also  be  noted.  These  features 
are  not  visible  in  the  optical  images.  Forest  regions  present 
texture  patterns  in  both  optical  and  radar  images  with  a  shift 
or  reversal  of  intensity.  Urban  structures  such  as  buildings 
do  not  appear  in  great  detail  in  the  radar  images  but  shadowing 
effects  are  apparent  in  the  radar  images.  The  predominant  simi¬ 
lar  features  in  both  optical  and  radar  images  are  shapes  of  cer¬ 
tain  features  such  as  roads,  forest  boundaries  and  land-water 
boundaries.  Furthermore,  these  shapes  remain  recognizable 


even  after  some  degree  of  degradations  in  resolution,  intensity, 
geometry  or  weather.  Thus,  it  may  be  expected  that  structural 
features  will  be  important  in  map  matching. 
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!•  i  H»i  rf*  2.  Typical  radar  iniaeo  cons  idc  red  in  this  study 


1.1  PROBLEM  STATEMENT  AND  BACKGROUND  OF  PREVIOUS  WORK 

Let  two  images,  S  the  search,  and  W  the  window,  be  as 
shown  in  Figure  3.  S  is  taken  as  an  MxN  array  of  digital 
picture  elements  which  assume  one  of  K  gray  levels: 

0  <  S  (m,  n)  <  K  -  1 
1  <  m  <  M 
1  <  n  <  N . 

W  is  taken  as  a  JxK  array  of  digital  picture  elements 
(with  J  <  M  and  K  <  N)  having  one  of  Q  gray  levels: 

0  <  W(j,k)  <  Q  -  1 
1  <  j  <  J 
1  <  k  <  K. 

When  W  is  superimposed  on  S,  the  position  of  the  upper  left  hand 
corner  of  W  is  used  as  a  reference  point  in  identifying  the 
position  of  W  in  the  search  area  coordinates.  Therefore,  when 
W  is  located  at  (u,v) ,  a  picture  element  S(i,j)  in  S  is  super¬ 
imposed  by  a  picture  element  W(i-u,j-v)  in  W. 

In  the  physical  world,  S  may  represent  a  high-resolution 
optical  picture  taken  by  an  air-borne  camera  at  high  altitude. 
This  is  as  depicted  in  Figure  4.  W  may  represent  an  image  taken 
by: 

.  The  same  sensor  at  a  different  time  and  look  angle 

.  A  difference  sensor  such  as  a  side-looking  radar  at 


a  lower  altitude 


FIGURE  3.  SEARCH  SPACE 


The  basic  problem  is  that  given  an  object  terrain  as 
represented  by  W,  we  wish  to  determine  by  the  use  of  a  digital 


computer  whether  the  same  object  terrain  appears  in  S. 

A  number  of  researches  have  been  made  to  develop  methods 
of  detecting  local  similarity  and  perform  image  matching.  Some 
of  the  more  promising  ones  are: 

.  Basic  correlator 
.  Statistical  correlator 
.  Sequential  correlator 
.  Sequential  template  matching 

In  the  following  sections,  each  of  these  methods  is  briefly 
described  and  their  relative  performances  are  analyzed. 

1.1.1  Basic  Correlator 

The  basic  correlation  is  a  method  used  to  form  a  correlation 
measure  between  two  picture  functions  and  to  determine  the 
location  of  the  maximum  correlation  [l],  [2],  In  applying  this 
technique,  the  correlation  measure  R(u,v)  at  the  reference 
location  (u,v)  is  defined  as: 


K 

E 


r>/  \  j  =  1  i  =  1 

R(u,  v)  =  J 


S  S{i,  j)  W(i-u,  j-v) 


(7” K  1  ;  n  r  K  3  2  ^ 

|  S  E  S  (i,  j)  I  !  2  E  W  (i-u,  j-v)lp 
<-j  =  l  i  =  l  -  -j  =  l  i  =  l 


(1) 


For  a  given  window  W, 


the  term 


ll  W2 (i-u, j-v) 

ji 


in  Equation  (1)  is 
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a  constant  for  any  reference  point  u,v.  Therefore  Equation  (1) 
can  be  reduced  to: 


R(u,  v) 


K  J 

^  ^  jS(  i>  j)  W( i-u,  j  -v) 


(2) 


To  determine  the  location  of  maximum  correlation,  R(u,v) 
must  be  computed  at  each  location  u,v;  1  <  U  <  (M-J+l) ,  1  <  V  < 
(N-K-l) .  This  is  because  no  decision  can  be  made  until  the 
correlation  array  R(u,v)  is  computed  for  all  u,v.  The  perfor¬ 
mances  of  this  correlator  can  be  described  as  follows: 

(1)  This  method  is  relatively  sensitive  to  image  noise  [3], 

In  the  presence  of  image  noise,  the  correlation  function  pro¬ 
duces  a  relatively  broad  peak,  thus  making  a  selection  of  a 
correlation  peak  difficult. 

(2)  A  great  amount  of  computation  must  be  performed  since 
the  window  and  search  areas  are  usually  large  in  an  actual 
photograph.  With  this  technique  no  decision  can  be  made  until 
the  correlation  array  R(u,v)  is  computed  for  all  u,v. 

(3)  With  the  exception  to  the  matching  of  every  simple 
pictures  this  method  and  the  variations  of  this  method  [4],  [5] 
does  not  provide  satisfactory  performance  in  image  matching. 

1.1.2  Statistical  Correlation 

To  overcome  some  of  the  difficulties  mentioned  in  the 
previous  method,  the  statistical  knowledge  of  the  spatial  relation 
of  picture  elements  within  each  image  were  used  in  this  statistical 
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correlation  [6].  The  statistical  correlation  measure  Rs(u,v)  is: 


Rs(u,v)  = 


K  J 

S  E  F  (i,j)F  (i-u,j-v) 

J  =  1  1~1 _ ^ 

r  k  j  ; k  j  ; Zf  **• 
rs  Sf^jM  S  Sr‘[i,j) 

uj  =  l  i  =  l  8  J  Lj  =  l  1  =  1  w  J 


(3) 


where  F  (i, j)  and  F  (i , j )  are  obtained  by  spatially  convolving 
s  w 

the  images  S  ( i , j )  and  W(i,j)  with  spatial  filter  functions  Ds  ( i , j ) 
and  Dw(i, j) : 


Fs(i,j)  =  S(i, j)  ©  Ds(i,j)... 

(4) 

Fw(i,j)  =  W(i,j)  ©  Dw(i,j)... 

(5) 

The  spatial  filter  functions  D  (i,j)  and  D  ( i , j )  are  chosen 

s  w 

to  maximize  the  correlation  peak.  The  first  step  in  the  spatial 
filter  design  is  to  decorrelate  or  whiten  the  images  as  follows: 


A  =  [Hs]-1S  (6) 

B  =  [H  ]_1W  (7) 

w 

S  and  W  are  column  vector  representations  of  the  images  S(i,j) 

and  W(i, j)  obtained  by  column  scanning  the  images.  H  and  H  are 

s  w 

obtained  by  a  factorization  of  the  image  covariance  matrices 


k  =  hh; 
s  s  s 


Kw  -  HwHw 


(8) 

(9) 


H  and  H  may  be  formulated  in  terms  of  the  eigenvector  and 
s  w 

eigenvalues  of  K  and  K  as  follows: 
s  w 


T*  4  4 

k=eae=(ea  me  a  r 

s  s  ‘s  s  '  s'V  '  s‘V 


=  h.h.t 


(10) 
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K  =  E  A  ET  =  (E  A  )**  (E  A  )  ** 
W  WWW  w  w  w  w 


T 

=  H  H 
w  w 


The  correlation  operation  is  performed  on  the  whitened 
vectors  A  and  B  to  yield  the  statistical  correlation  measure. 


Rg(u,v)  = 


ATB 


(ata)Js(btb)J5 


or 


where 


RSK  v)  = 


T  “1  T 

(K  v  s  W 


{[(K1)"  S]T[(KT)'1S]WTfj5 


(KT)_1  =  (HsH^)-1, 


If  the  image  elements  are  assumed  to  be  samples  of  Markov 
process  then: 


n 

1 

-pS'1  0 

0 

1  o 

T  -1  -1  ,  2 

-P2'1 

2-1  -1 
(1+P  )E  -P2 

0 

. .  .  0 

(K  )  =  K  =  1/(1  -p  ) 

0 

-1  2-1 
-pE  (1+p  )Z 

-pS'1 

.  .  .  0 

0 

«  *  4 

1 

(15) 


where 


(ID 


(12) 


(13) 


(14) 


r  =  correlation  between  adjacent  image  elements 


E_1=  l/(l-p2) 


-p  0  0 


-p  (1+p  )  -P  o 


-p  0  0 


•  .  T  —  1 

Multiplying  S  by  the  (K  )  is  equivalent  to  convolving  the 
image  S(i,j)  with  the  spatial  filter  function  D(i,j). 


D(i.j)  = 


-p(l+p2) 


-p(l+p  ) 
(1+p2)2 
.-P(l+P2) 


-p(l+p2) 


and  Equation  (4)  and  (5)  become 


Fs(i/j)  =  S(i,j)  8  D  (i , j ) 


Fw(i, j )  =  W(i, j )  @  D (i, j ) 


Performance  of  this  correlator  can  be  described  as  follows: 

(1)  Experiments  performed  on  selected  images  indicated  that 
the  statistical  method  does  provide  better  performance  in  terms  of 
providing  a  sharper  peak  at  the  location  of  image  matching.  In 
order  for  the  method  to  work  well,  prior  knowledge  of  the  picture 
statistics  is  required.  In  an  actual  situation,  this  information  is 


either  not  available  or  extensive  computations  are  required  so  that 
the  correlator  can  be  re-designed  to  tailor  to  the  input  data. 

(2)  It  is  anticipated  that  performance  of  this  correlator 
would  be  worse  than  that  of  the  basic  correlator  if  the  statistics 
of  the  input  data  differ  from  statistics  used  in  the  design  of  the 
correlator . 

(3)  Because  of  extra  convolution  steps  needed  in  this  method, 
the  amount  of  computations  needed  by  this  method  is  even  more  than 
that  needed  for  the  basic  correlator.  It  remains  to  show  that  the 
improvement  in  performance  is  worth  the  extra  cost  in  computations. 

1.1.3  Sequential  Correlator 

A  common  criticism  of  both  the  basic  and  statistical  correlators 
is  the  great  amount  of  computations  that  must  be  performed.  A 
method  of  sequential  correlation  has  been  proposed  [7]  to  reduce 
the  computation  time.  The  basic  form  of  this  algorithm  is  simple. 

An  error  function  is  defined  as  follows: 

K  J 

5{u,  v)  =  E  E  |S(L,  j)  -  W(i-u,  j-v)|  ..»«  (20) 

j  =  l  i  =  1 

Instead  of  testing  each  of  the  elements  in  the  window  area, 
elements  of  the  area  are  selected  at  random.  The  error  is  accumu¬ 
lated  for  as. each  of  the  elements  is  compared.  If  the  error  exceeds 
a  predetermined  threshold  value  before  all  the  elements  in  the 
window  area  are  tested,  the  test  is  considered  failed  for  the  window 
(u , v )  and  a  new  window  is  tested.  The  test  procedure  is  depicted 
in  Figure  5.  Curves  A,  B,  and  C  depict  the  cumulative  errors  for 
three  different  reference  points.  A  and  B  accumulate  errors  rapidly 


and  the  tests  terminate  early.  Curve  C,  however,  accumulates 
error  more  slowly.  It  is  therefore,  much  more  likely  to  be  a 
candidate  as  a  matching  point.  Theoretical  analyses  and  simulation 
tests  [7]  indicated  that  a  saving  of  computation  time  of  at  least 
two  orders  of  magnitude  is  possible.  As  with  the  other  two 
correlation  methods  mentioned  earlier,  only  very  simply  structured 
images  were  processed  successfully  by  this  new  approach. 

1.1.4  Sequential  Template  Matching 

Extending  the  basic  concept  described  in  1.2.3  a  more  compli¬ 
cated  sequential  detecting  method  was  proposed  [8],  Local 
similarity  between  a  set  of  templates  is  matched  to  a  given  image. 
Instead  of  matching  each  template  of  a  set  to  an  image  at  every 
location,  the  templates  are  partitioned  and  a  representative 
template  is  defined  for  each  of  the  partitions.  Several  levels 
of  partitions  are  defined.  Elimination  of  mismatching  locations 
and  termination  of  computation  can  take  place  at  each  level  of 
detection.  Each  level  of  testing  is  over  a  more  restrictive 
subset  of  template  class  than  the  previous  level.  Matching  process 
terminates  when  the  accumulative  template  matching  error  exceeds 
a  threshold  level.  A  location  which  has  gone  through  successive 
levels  of  matching  without  rejection  is  declared  a  likely  candidate. 

The  performance  of  the  sequential  template  matching  method  can  be 
described  as  follows: 

1.  Computation  time  is  reduced  due  to  the  sequential  method 
of  testing.  In  matching  a  real  image,  several  hundreds  of  tem¬ 
plates  are  needed,  thus  making  the  task  of  template  partition  difficult. 

2.  In  a  sequential  testing,  the  ordering  of  the  features  is 
important.  Essential  features  (such  as  roads,  rivers,  etc.)  must 


appear  in  both  images  in  order  to  be  considered  as  a  match. 
Methods  of  ordering  features  have  not  yet  been  developed. 


Tgui—MniNninwg 


18 


1.2  OUTLINE  OF  REPORT 

Two  general  approaches  to  map  matching  may  be  differentiated. 
The  first  approach  may  be  called  a  pictorial  method  and  generally 
involves  a  transformation  of  one  image  into  a  new  image  which 
corresponds  to  the  reference  image.  The  second  approach  may  be 
called  a  feature  method  of  the  emphasis  is  placed  on  the  second 
more  robust  method. 

The  pictorial  method  is  developed  in  Section  II.  Transfor¬ 
mations  for  matching  map0  from  different  sensors  and  geometries 
are  analyzed.  Examples  of  perspective  transformations  are  given 
for  optical,  radiometric,  and  radar  images.  A  sensor  transfor¬ 
mation  based  upon  a  statistical  matching  of  corresponding  image 
points  is  also  presented. 

The  use  of  edge  features  for  map  matching  is  described  in 
Section  III.  A  new  orthogonal  decomposition  in  terms  of  point, 
edge,  line,  and  other  basis  vectors  is  presented.  The  advantage 
of  this  representation  in  a  local  region.  Map  matching  using 
edge  features  also  permits  high  speed  logical  similarity  mea¬ 
surement.  The  number  of  computations  using  this  technique  is 
especially  attractive  for  real  time  implementations.  The  use  of 
edge  features  for  map  matching  is  also  demonstrated.  In  the  ab¬ 
sence  of  geometric  and  sensor  distortions,  the  method  works  well 
as  indicated  by  the  results  of  the  experiments  on  matching  an  edge 
image  to  a  noise-corrupted  version  of  "the  same -image.  For  optical 
to  radar  matching,  the  performance  is  somewhat  degraded  but  can 
be  improved  using  geometric  and  sensor  intensity  corrections  as 
well  as  data  pre-processing* 
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The  theory  of  invariants  is  reviewed  in  Section.  The 
study  of  algebraic  invariants  is  an  important  field  of  mathe¬ 
matics  and  several  results  of  this  theory  are  summarized. 

The  concept  of  perceptual  invariants  has  also  been  studied 
for  the  visual  system.  Several  invariant  properties  of  the 
visual  system  are  presented  and  a  simple  device  which  produces 
invariant  measurements  is  described.  The  use  of  spatial  mo¬ 
ments  for  invariant  measurements  is  described  in  Section  4.3. 

A  powerful  theorem  which  relates  moment  invariants  to  alge- 
!  braic  invariants  for  continuous  functions  i.?.  stated.  New 

results  indicating  a  degree  of  variance  in  the  invariant 
measurements  for  discrete  functions  are  presented.  Examples 
are  given  which  also  illustrate  that  the  magnitude  of  the 
variance  can  be  controlled  by  careful  processing. 

Sequential  techniques  for  searching  for  matching  image 
locations  are  presented  in  Section  V.  These  methods  promise 
logarithmic  efficiency  over  conventional  correlation  techni¬ 
ques.  Zoom  techniques  are  an  important  special  case  of  the 
sequential  techniques  and  are  considered  in  detail.  The 
effects  of  reqion  size  and  resolution  limitations  are  con¬ 
sidered  and  demonstrated  by  several  examples.  Finally,  recom¬ 
mendations  for  future  work  in  map  matching  are  given  in 


Section  VI. 
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SECTION  II.  TRANSFORMATIONS  FOR  MAP  MATCHING 

2.1  GEOMETRIC  TRANSFORMATIONS 

In  this  section,  methods  of  registration  of  the  optical 
and  radar  images  are  described.  Figure  6  shows  the  general 
geometry  of  the  optical  and  radar  imaging  systems  together 
with  an  object  terrain  of  interest.  In  this  figure,  the 
optical  system  is  assumed  pointing  at  the  nadir  and  the  radar 
system  is  scanning  the  same  object  terrain  at  a  slant  range. 
Figure  7 (a)  depicts  a  rectangular  array  of  digitized  pixels 
of  the  optical  image.  Figure  7  (b)  represents  the  rectili- 
nearly  digitized  radar  image  plotted  in  the  optical  image 
coordinates.  Point  (x,y)  in  the  optical  image  and  point  (x, 
y)  in  the  radar  image  represent  an  image  of  the  same  object 
point  in  the  terrain  of  interest.  The  basic  transformation 
consists  of  deforming  the  radar  image  so  that  its  features 
correspond  as  closely  as  possible  to  those  in  the  optical 
image . 

Three  methods  of  transformations  were  studied.  These 
are:  •  Perspective  Transformation 

•  Polynomial  Estimate 

•  Interactive  Nonlinear  Transformation 

2.1.1  Perspective  Transformation 

The  general  geometry  of  the  optical  imaging  system  is 
shown  in  Figure  8.  In  this  figure,  the  object  (ground) 
plane  is  assumed  to  be  parallel  to  the  image  plane.  The 


optical  axes  is  perpendicular  to  both  planes  and  the  en¬ 
trance  pupil  of  the  camera  is  located  at  the  origin  of  a 
rectangular  coordinate  system.  From  the  lens  law  of  geo¬ 


metrical  optics,  the  object  distance  u3 ' 
are  related  by: 


and  image  distance 

(1) 


where  f  is  the  optical  system  focal  length.  Since  the 
object  distance  is  large  compared  to  f,  the  image  plane  is 
the  fecal  plane  and  x^  =  f.  Any  object  point  u'  =  (u^  ' , 
u3 ' /  u3 ' )  is  imaged  into  a  point  x  =  (x^,  x^,  x^)  in  the 
focal  planes.  The  components  of  x  can  be  computed  by  the 
following  equations: 


XlSf 


u  • 
1 


X 


2 


i 


(2) 


(3) 


(4) 


The  geometrv  of  the  radar  imaging  system  is  shown  in 
Figure  9.  In  this  geometry  the  radar  image  plane  is 
assumed  to  be  rotated  0,  0  and  4)  in  the  pitch,  roll  and  yaw 
axes.  Therefore,  the  orientation  of  the  radar  system  is 
specified  by  the  angles  of  rotation  0,  9  and  <p  .  An  object 
with  a  ground  coordinates  u  =  (u^,  u^,  u^)  in  the  radar 


imaging  system  can  be  related  to  the  old  coordinates  u' 


(u^,  u^,  u^)  in  the  optical  imaging  system  by  the  matrix 
equation 

u*  =  Mu  (5) 
where  M  is  the  product  of  three  rotational  transformation 
matrices : 
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Successive  rotations  in  a  pointing  sequence  can  be 
defined  by  M  which  is  the  product  of  thr  three  matrices 
(6),  (7),  and  (8).  A  list  of  these  product  matrices  for 

all  six  permutations  of  1-,  2-,  and  3-  axis  are: 


(cos  0  cos  0)  (cos  0  sin  0  sin  cp  +  sin  0  cos  Cf)  (cos  8  sin  0  cos  cp  +  sin  0sincp) 

(sin  0cos  0)  (cos  0cos  cp  +  sin  0s  in  0s  in  cp)  (sin  0  sin  0  cos  cp  -  cos  0sincp) 

(-sin  0)  (cos  0  sin  cp)  (cos  0  cos  cp) 

(cos  0  cos  0)  (sin  0  sin  cp  -  sin  0  cos  0  cos  cp)  (sin  0  cos  0  sincp+  sin  0  cos  cfi 

(sin  0)  (cos  0  cos  cp)  (-cos  0  sin  cp) 

(-cos  Osin  0)  (sin  8  sin  0  cos  cp  +  cos  0  sin  cp)  (cos  0  cos  cp  -sin  Osin  0  sin  cp) 
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(cos  P  cos  0  -sin  P  sin  0  sin  cp)  (-sin  0  cos  cp)  (sin  {3  cos  0  +cos  0  sin  0  sincp) 
(sin  P  cos  0  sin  cp  +  cos  p  sin  0)  (cos  0  cos  cp)  (sin  P  sin  0  -cos  P  cos  0  sincp) 
(-sin  p  cos  cp)  (sin  0)  (cos  0  cos  cp) 


m23i  = 


(cos  0  COS  P) 

(sin  0  cos  p  cos  cp  +  sin  P  sincp) 
(sin  0  cos  P  sincp  -  sin  p  cos  cp) 


(-sin  0)  (cos  0  sin  p) 

(cos  0  cos  cp)  (sin  0  sin  p  cos  cp  -  cos  p  sin  cp) 
(cos  0  sincp)  (cos  p  cos  cp  +  sin  0  sin  P  sincp  ) 


m 


312' 


m 
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(cos  0  cos  P  +  sin  0  sin  P  sin  cp) 
(sin  0  cos  cp) 

(sin  0  cos  P  sincp  -cos  0  sin  p) 


(cos  0  sin  P  sincp  +  sin  0sin  0)  (sin  P  cos  cp) 
(cos  0  cos  cp)  (-sin  cp)(-sincp) 

(cos  0  cos  P  sincp  +  sin  0  sin  P)  (cos  Pcoscp) 


(cos  0  cos  P) 

(sin  0  cos  cp  +  cos  0  sin  P  sincp) 
(sin  0  sincp  -  cos  0  sin  0  cos  Cf) 


(-sin  0  cos  0)  (sinp) 

(cos  0  cos  cp  -  sin  0  sin  P  sincp )  (-cos  p  sincp) 
(sin  0  sin  0  cos  cp  +  cos  0  sincp)  (cos  p  cos  cp)  _ 


The  matrices  (6),  (7)  and  (8)  and  their  products  M  (in 

any  order) are  orthonormal,  nonsingular,  and  have  unity  deter¬ 
minates.  These  properties  imply  that  the  transform  in  (5) 
can  be  inverted  as: 

u  =  M'1  u1  (9) 

Substituting  (5)  into  (2),  (3)  and  (4)  yields: 


u. 


u3)  =  f 


milUl  +  m12U2*ml3U3 
m31Ul  +  m32U2  +  m33  U3 


(10) 


.-1 


X2  =  C2  (U1  ’  u2  ’  U3}  =  f 


m21Ul  +  m22  U2  +  m23U3 


m31  U1  +  m32  U2  +  m33  U3 


(ID 


& 


x3  =  C3  (u-^  u2,  u3)  =  f  (12) 

where  rru  j  are  the  elements  of  the  product  matrix  M. 

Let  the  radar  imaging  system  be  located  at  an  altitude 
h.  Substituting  h  for  u3  in  Equations  (10,  (11)  and  (12), 

we  have: 


*1  =  CI  <U1  •  u2>  *  mXl  ul  4  ml2  “2  j  ml3h 

m31  U1  +  m32  U3  +  m33h 


*2  =  c\  (ul'  U2>  =  m21  U1  +  m22  U2  +  ">23h 

m31  "l  +  m32  "2  +  m33h 


*3 =  c2  (ur  U2>  =  f 


In  the  case  that  the  side  looking  radar  is  tilted  at  an 
angle  3  and  «J>  =  e  =  0,  equations  (13)  and  (14)  become: 


x  =  f 


-Uj  cos  p  +  h  sin  P 
Uj  sin  p  +  h  cos  P 


c2  ~  Uj  sin  p  +  h  cos  P 
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Given  the  grid  points  (u^,  u^)  in  the  object  plane 
coordinates.  Equations  (13)  -  (17)  can  be  used  to  compute  the 
corresponding  points  (x^,  X2)  in  the  radar  imaging  system 
coordinates.  in  a  quantized  system,  the  computed  coordi¬ 
nate  points  do  not  generally  fall  on  the  grid  points  in  the 
radar  image.  So  it  is  necessary  to  adapt  rules  for  picking 
the  appropriate  element  out  of  the  radar  image.  The  easiest 
method  is  to  simply  pick  out  the  value  of  the  nearest  neigh¬ 
bor  to  the  computed  grid  point.  This  method  works  well  when 
the  geometric  distortion  is  not  severe  and  when  the  x  and  u 
grids  are  about  the  same  size.  \  second  method  is  to  take 
the  average  (possibly  weighted  average)  of  the  neighbor  ele¬ 
ments.  A  third  and  more  complicated  method  is  to  use  inter¬ 
polating  polynomials  as  a  means  of  computing  values  of  the 
elements  between  points  of  known  values. 

2.1.2  Polynomial  Estimate 

Let  the  optical  image  intensity  be  given  by  f(u^,  Uj) 
and  the  radar  image  intensity  be  given  by  f (x^,  ^2^  ‘  A 
technique  can  be  used  to  map  the  radar  image  into  the 
optical  image  using  the  following  equations: 


=  u2) 

(18) 

=  g2<V  u) 

(19) 

g^  and  g2  can  be  approximated  by  a  polynomial  in  u^  and  U2 
of  degree  N  as: 

N  N-i  i  j 

*1  -  ci  <v  V  ’  EE  hu  ”i°z 

i=o  j=o 


(20) 


(21) 


N  N-i 

x  =  C  (»,*  u2^  “ 

L  i=o  3=0 


i 


where  k. . .  and  k„. .  are  the  constant  polynomial  coefficients. 
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For  most  practical  cases,  a  second  degree  (N  =  2)  approxi¬ 
mation  is  adequate.  This  has  been  verified  by  experi- 
ments[9]  (on  data  from  side  look  radar  iamges)  in  which  it 
was  found  that  the  decrease  in  rms  error  as  a  function  of 
increasing  polynomial  terms  is  not  significant  for  N  greater 
than  2 . 

The  coefficients  of  Equations  (20)  and  (21)  can  be 
computed  by  fitting  the  two  dimensional  functions  to  a  set 
of  f(u1#  u2)  and  f(x1#  x2)  values.  The  linear  -  least  - 
squares  estimate  procedure  can  be  used  to  express  C1(u1,  u2) 
surfaces  by  polynomials  whose  squared  distances  from  the 
true  surface  is  a  minimum.  To  obtain  the  coefficients,  pro¬ 
minent  features  that  appeared  in  both  the  optical  and  radar 
images  are  selected.  These  prominent  features  can  be: 

•  End  points  of  a  long  edge 

•  Intersection  of  two  lines 

•  Corners  of  detected  squares 

•  Distinguishable  points  on  the  boundary  of 
texture  region 

For  a  second  degree  (N  =  2)  approximation,  at  least  six 


pairs  of  conjugate  points  are  needed.  The  values  of  these 
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points  can  be  arranged  in  the  following  manner: 
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2  2 


*2  U1U2 
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- 

(22) 


Where  R  •>  6. 

A  similar  equation  exists  for  x2  and  k2ij'  BotVl  of  these 
may  be  written  in  matrix  notation  as: 


xx  =  U  K1 

(23) 

x2  =  0  K2 

(24) 

Using  the  linear 

-  least  -  squares  estimate 

theory. 

the  best  estimate  for 

and  K2  are  given  by  the 

pseudo 

inverse  solution: 

T  -1  T 

(25) 

K1  = 

:  (U  U)  U  X1 

T  -1  T 

K2  = 

:  (U  U)  U  X2 

(26) 

2.1.3  Interactive  Nonlinear  Transformation 

In  this  section,  a  method  of  removing  spatial  distor¬ 
tion  in  the  radar  imaqe  with  respect  to  the  optical  image 
is  described.  This  method,  which  was  first  described  by 
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Ulstad  [10]  is  currently  being  studied  as  a  possible  method 
of  removing  geometric  distortion  by  a  nonlinear  spatial 
transformation. 

The  spatial  registration  system  is  shown  in  Figure  10. 
The  two  sets  of  image  data  are  sequentially  available  as 
columns  of  matrices  X  and  U  each  with  L  rows.  The  L-element 
column  vectors  of  X  and  U  are  processed  by  the  system  in 
order  to  produce  a  third  L-element  column  vector  R.  The  L 
elements  of  the  comumn  of  U  map  into  L  corresponding  points 
in  X  based  on  the  following  equation. 

u(i,k)  =  x  ^i  +  6  (i,k)  ,  k  +  a  (i,k)j  (27) 

The  quantities  B(i,k)  and  a  (i,k)  are  real  scalar  cor¬ 
rection  factors  are  functions  of  the  indices  N  and  k. 

As  shown  in  Figure  11,  use  is  made  of  N  submatrices 

defined  in  U.  Submatrix  U  consists  of  2  p  rows  labeled 

n 

sn  -  p  through  sn  +  p  where  1  <_  n  <_  N,  p  is  an  integer  and 


Let  the  N  index  points  (a  b  )  in  Matrix  X  be  such 

n  n 

that: 

U (sn,k)  =  X(an,  bn) ,  (29) 

The  points  (an,  t>n)  in  X  correspond  to  the  middle  row 

of  the  nth  submatrix  U  at  the  kth  column.  The  points  will 

be  used  as  the  corner  points  of  a  piecewise  linear  synthetic 

scan  line  in  the  transformation.  In  determining  (a  ,  b  ) 

n  n 

during  the  processing  of  the  kth  column,  correlations  of 
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FIGURE  10.  SPATIAL  REGISTRATION 


radar  image  X 


compute 

corrections  in 
corner  points  of 
synthetic  scan 


optical  image  U 


compute 

synthetic 

scan 


compute 

correlations 

-transformed  radar  image,  R 


FIGURE  11.  INTERACTIVE  TRANSFORM  SYSTEM 
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R(i,k)  of  matrix  R  and  U(i,k)  of  matrix  U  are  made  for  i, 
in  ranging  from  sn-p  to  sn+p.  The  results  of  these  corre¬ 
lations  is  then  used  to  compute  a  new  corner  ooint  (a  + 

n 

Aa  ,  b  +  Ab  ) .  The  corner  point  with  the  updated  value  can 
n  n  n 

then  be  used  to  produce  the  (k+l)th  column  of  the  R  matrix. 
This  process  continues  until  an  entire  matrix  R  is  produced. 

2.1.4  Results 

The  polynomial  method  is  preferred  over  the  perspective 
transformation.  Perspective  transformation  requires  tracking 
and  ephemeris  data  to  determine  the  positions  of  the  sensors 
at  the  time  of  exposure.  These  data  are  often  not  avail¬ 
able.  In  the  polynomial  method,  data  needed  for  the  trans¬ 
formation  are  entirely  contained  in  the  images. 

Figures  12,  13  and  14  show  the  results  of  geometric 
distortion  corrections  by  the  method  discussed  in  Section 
2.1.2.  Figure  12  shows  the  various  possible  degrees  of 
corrections.  Figure  13  shows  a  radiometric  and  an  optical 
picture.  Geometric  transformation  was  performed  on  the 
optical  picture  in  order  to  match  the  radiometric  picture. 
Corrections  were  concentrated  on  matching  the  three  straight 
roads.  A  total  of  eight  reference  points  from  each  of  the 
two  original  pictures  was  used  for  the  transformation. 

Figure  14  shows  two  pictures  taken  by  two  different  sensors, 
an  optical  camera  and  a  radar,  at  different  look  angles. 

The  transformation  is  made  on  the  optical  image  so  that  it 
can  be  geometrically  in  registration  with  the  radar  image. 


Radar 


Mi  c rad 


Transformed  Mi c rad  Transformed  Micrad 


Figure  12.  Polynomial  Geometric  Corrections 
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For  this  transformation,  a  total  of  twelve  reference  points 
from  each  of  the  two  original  images  is  used  to  form  based 
for  a  pseudo -inverse  transformation.  As  shown  in  Figure  14 
(d) ,  a  nearly  perfect  registration  was  obtained. 


t 


2.2  Sensor  Transformation 
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Because  of  the  differences  in  operating  characteristics 
of  the  two  sensors,  images  of  the  same  object  taken  by  the 
optical  system  and  the  radar  system  will  have  different  in¬ 
tensity  profiles.  The  most  prominent  difference  is  that  the 
optical  image  is  positive  while  the  radar  image  is  negative. 
In  this  section,  methods  of  intensity  transformation  are 
described.  These  methods  consist  of  transforming  the  radar 
image  so  that  its  intensity  profile  matches  as  closely  as 
possible  to  that  of  the  optical  image.  The  two  methods 
studied  are: 

•  Intensity  Equalization  using  Karhunen-Loeve  Trans¬ 
form 

•  Equalization  using  Intensity  Averaging 
2.2.1  Intensity  Reversal 

As  the  first  step  in  the  intensity  transformation,  the 
negative  radar  image  is  transformed  into  a  positive  image. 
This  is  done  by  replacing  the  amplitude  of  each  picture 
element  e^j  by  its  complementing  value  e^.  For  an  image 
which  has  been  digitized  with  n  bits: 


where  2n  is  the  number  of  quantization  level.  Figure  14(a) 
shows  the  original  radar  image.  An  intensity  reversal  is 
performed  on  this  image  by  applying  Equation (30) .  The  result 
is  shown  in  Figure  15(a).  It  can  be  seen  that  the  contrast 
of  the  image  is  low  and  intensity  modification  must  be  per¬ 
formed  so  that  its  intensity  profile  matches  as  closely  as 


possible  to  that  of  the  optical  image  shown  in  Figure  15(b). 


I 


(a)  Radar  Image  (b)  Optical  Image 

Intensity  Reversed 


(c)  Intensity  Transformation  (d)  Intensity  Transformati 

(Karhunen-Loeve  Transform)  (Averaging) 


Figure  15.  Intensity  Transformation 
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2.2.2  Intensity  Equalization  Using  Karhunen-Loeve  Transform 
A  two-dimensional  image  intensity  function  f(x^,  x2) 
can  be  constructed  as  follows: 

f  ,  x2)  =  n  (31) 


where 

x^  =  Intensity  level  of  the  picture  element  in  the 
radar  image , 

0  £xx  <_2n  -1 

x2  =  Intensity  level  of  the  picture  element  in  the 
optical  image, 

0  <x2  <2n  -1 

n  =  The  numbers  of  picture  elements  which  have  inten¬ 
sity  level  of  x^  in  the  radar  image  and  the  inten¬ 
sity  level  of  x2  in  the  optical  image. 

A  covariance  matrix  E  is  then  formed  to  indicate  the 
relative  intensity  correlation  between  the  two  images  as 
follows: 

£t  =  E[(X1-Mxl)(X2-Mx2)T]  <32) 

where 


=  Expected  values  of 

Mx2  =  Expected  values  of  X2 

A  Karhunen-Loeve  transform  can  then  be  nerformed  to 
obtain  a  new  set  of  coordinates  and  as  follows: 


Zx  =<P  ZX<P  =  A 


(33) 


* 
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where 


A  -  Diagonal  matrix  with  the  eigenvalues  X.  of  2  as  the  diagonal 
elements 


cp-  Matrix  whose  columns  are  the  ordered  eigenvectors  of 


«P=  [«Pj.  «P2J  for  *  X2 


(34) 


2 

x 


The  transformed  coordinates  <t>  ^  and  7  are  shown  in 
Figure  16.  Based  on  and  <f>  2  intensity  corrections  can  now 
be  made  by  scaling  the  intensity  profile  of  the  radar  image 
such  that  the  resulting  coordinates  and  <}>  ^  are  orthogonal 
and  the  angle  between  and  the  f  ^  -  axis  is  approximately 
45  degrees. 

Figure  17  shows  a  three-dimensional  plot  of  the  resul¬ 
ting  intensity  transformation.  It  can  be  seen  from  this 
figure  that  both  images  have  approximately  the  same  numbers 
of  picture  elements  at  each  intensity  level.  Figure  15(c) 
shows  the  resulting  transformed  radar  image.  In  this 
figure,  considerable  details  which  are  not  visible  in  the 
original  radar  image  have  been  made  visible  by  the  inten¬ 
sity  transformation. 


2.2.3  Intensity  Equation  Equalization  by  Intensity  Averaging 

v 

The  intensity  transformation  method  described  in  Section 
2.2.2  can  be  modified  such  that  the  average  intensity  of  the 
transformed  radar  image  is  made  approximately  equal  to  the 
average  intensity  of  the  optical  image.  The  resulting  trans¬ 
formation  based  on  this  method  is  shown  in  Figure  15(d).  In 
this  method,  a  compromise  is  made  to  match  the  overall  in¬ 
tensity  of  the  two  images  whereas  image  transformation  using 


f  (Optical  ) . 

*•  *  \ 


FIGURE  16.  INTENSITY  COORDINATE  TRANSFORMATION 
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the  method  described  in  Section  2.2.2  products  a  better 
match  in  tracking  the  variations  of  local  intensity.  Matching 
the  overall  intensity  is  most  important  whenever  a  non- 
normalized  similarity  measure  is  used. 


2.3  Fourier  Transform 


The  most  salient  of  the  optical  systems  was  developed 
by  Lendaris  and  Stanley  for  terrain  classification.  The 
sensor  or  feature  extractor  measured  aspects  of  the  optical 
Fourier  transform  in  specified  geometric  patterns.  The 
initial  system  utilized  coherent  optical  techniques  to  ex¬ 
tract  measurements  (features)  of  the  Fraunhofer  diffraction 
pattern  which  were  classified  via  an  interactive  nonstatis- 
tical  classification  method.  In  the  present  application  the 
image  f(x,y)  is  transformed  optically  into  its  Fourier  space 
F(V  Vy)  via  the  well-known  property  that  coherent  light 
passed  through  an  image  f(x,  y)  and  then  through  a  thin  con¬ 
vex  lens  produces  the  diffraction  pattern  in  the  back  focal 
plane.  The  coordinates  (u,  v)  of  this  transform  plane  are 
directly  related  to  spatial  frequencies  and  of  the  input 
image.  This  is  shown  in  Eqs .  (35)  and  (36). 


oo  oo 

F(u,  v)  =  C 

-oo  -oo 


e 


-  2  tt  £  (ux  +  vy) 

Xf 


f(x,  y  )dxdy 


(35) 


= 


u 


“XT’ 


v 

y 


V 

tt- 


(36) 


The  constant  X  is  the  wavelength  of  the  incident  light,  f  is 
the  focal  length  of  the  transform  lens,  and  C  is  the  ampli¬ 
tude  transmission  factor.  The  sensor  is  composed  of  32 
angular  wedges  which  sample  the  squared  modulus  of  the  dif¬ 
fraction  pattern  (power  spectrum  I(u,  v) )  and  produces 
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appropriate  signatures 


I(u,  v)  =  F(u,  v)  F*  (u,  v)  =  |F(u,  v)J2  ,  (37) 

where  ~  denotes  a  complex  function  and  *  denotes  conjugate. 
The  annular  ring  sample  signature  is  expressed  in  polar 
coordinates  as 


r. 

J 


I(  P  0)  pd  pd0 


(39) 


j  =  1,2,. ...  32 
2  2  1/2  -1 

where  p  =  (u  +  v  )  and  8  =  tan  v/u.  Each  r ^  signa¬ 
ture  value  represents  the  total  power  in  the  annular  region 
(Pj»  Pj  +  Ap ) .  The  wedge  signature  values  may  be  similarly 
described  by 


Pmax  0.  +  £  e. 

Wj  ~  J  f  J  J  I(p,  0)pdpdfi  (39) 

Pmin  '  ^ 

j  =  1 ,  ....  32 

in  which  each  w^  represents  the  total  power  in  an  angular 
frequency  band  from  (9^,  8^  +  A9 )  . 

The  focal  length  f  of  the  transform  lens  is  adjusted  so 
that  the  outer  ring  of  the  detector  images  8  lp/mm.  The 
annular  ring  signature  provides  a  rotationally  insensitive 


50 


means  of  detecting  one-  or  two-dimensional  spatial  perio¬ 
dicity  within  apertured  areas.  In  contrast,  the  wedge  signa¬ 
tures  are  insensitive  to  periodic  structure  but  sensitive  to 
direction.  both  frequency  signatures  subsequently  have 
yielded  valuable  features  for  object  detection. 

To  initially  evaluate  the  effectiveness  of  Fourier 
features  for  map  matching,  the  following  experiment  was 
performed.  Four  corresponding  regions  in  an  optical  and 
radar  scene  were  selected  as  illustrated  in  Figure  18.  The 
Fourier  transform  of  the  corresponding  regions  were  computed 
using  a  coherent  optical  system.  Comparative  photographs  of 
the  magnitudes  of  the  Fourier  transforms  of  the  regions  are 
shown  in  Figures  19  to  22.  Note  that  the  most  obvious  pat¬ 
tern  is  that  the  optical  image  contains  more  high  frequency 
energy  corresponding  to  greater  detail  in  the  optical  images. 
Also,  a  correlation  between  line  structures  in  the  images 
and  transforms  is  visible  in  several  of  the  patterns.  To  in- 
vestiqate  the  transforms  with  a  more  sensitive  detector, 
ring  and  wedge  measurements  were  made  using  the  RSI  detector. 
Graphs  of  these  results  are  shown  in  Figures  23  to  30.  Note 
that  these  measurements  indicate  a  correlation  in  both  ring 
and  wedge  measurements.  Thus,  more  quantitative  studies 
such  as  map  matching  with  these  measurements  are  required. 


(a)  Fourier  transform  of 
optical  image  of  the 
rural  field  region 


(b)  Fourier  transforn'i  of 
radar  image  of  rural 
field  region 


Fourier  transform  of  rural  field  region 


(a)  Fourier  transform  of  optical 
image  of  road  intersection 
region 


(b)  Fourier  transform  of  radar 
image  of  road  intersection 
region 


Fourier  transform  of  road  intersection  reg 
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(a)  Fourier  transform  of  optical 
image  of  curving  road  region 


♦ 


(b)  Fourier  transform  of  radar 
image  of  curving  road  region 


Figure  21. 


Fourier  transform  of  curving  road  region 


(I) )  Fourier  transform  of  radar 
image  of  stadium  region 


Figure  22. 


Fourier  transform  of  the  stodium  region 


FOURIER  TRANSFORM  SIGNATURE  OF  THE.  OPTICAL  IMAGE  OF  RURAL  FIELD  REGION 
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Figure  25. 

FOURIER  TRANSFORM  OF  OPTICAL  IMAGE  QF  ROAD  INTERSECTION  REGION 


(lp/mm):  0. 0  1. 21 


FOURIER  TRANSFORM  OF  RADAR  IMAGE  OF  ROAD  INTERSECTION 


(Ip /mm):  0.0  1. 21 


FOURIER  TRANSFORM  SIGNATURE  OF  THE  OPTICAL  IMAGE  OF  THE  CURVING  ROAD  REGION 


(lp/mm):  0.  0  1. 21 
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Figure  29. 

FOUIER  TRANSFORM  OF  OPTICAL  IMAGE  OF  STADIUM  REGION 


(lp./mm):  0. 0  1. 21 
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SECTION  III.  EDGE  EXTRACTION  AND  REGISTRATION 

3 . 1  Introduction 

An  important  technique  used  in  automatic  scene  analysis 
is  image  segmentation.  This  technique  involves  separating 
the  image  into  various  regions  corresponding  to  individual 
objects.  Assuming  that  these  regions  have  some  homogeneous 
characteristic ,  for  example,  luminance,  color,  texture,  etc., 
one  segmentation  technique  is  to  detect  sharp  transitions 
called  edges,  which  tend  to  outline  the  desired  boundaries. 
The  opposite  alternative  is  to  "grow"  regions  by  connecting 
small  adjacent  areas  of  similar  characteristics.  Of  interest 
here  is  the  detection  of  edges  that  separate  regions  of  dif¬ 
ferent  constant  luminances,  and  lines  which  can  be  regarded 
as  a  degenerate  pair  of  edges.  This  operation  requires  the 
examination  of  several  picture  elements  within  contiguous  or 
overlapping  sub-areas  of  the  image,  followed  by  a  decision  as 
to  whether  an  edge  or  the  line  segment  is  present  or  not 
within  each  sub-area. 

The  segments  can  be  characterized  by  variable  such  as 
amplitude,  orientation,  position  within  the  sub-area,  etc., 
and  possibly  a  measure  of  confidence. 

Upon  examination  of  the  whole  picture,  the  object  boun¬ 
daries  are  constructed  by  connecting  the  edge  and  line  ele¬ 
ments  detected  previously.  This  operation  can  be  directed  by 
simple  syntactic  rules,  for  example  connect  neighbor  edge 
elements  that  line  uo  approximately,  and  delete  isolated 
parallel  elements. 
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Some  of  the  difficulties  of  edge  detection  are  caused 
by  noise,  but  much  more  so  by  the  fact  that  visually  dis¬ 
tinct  edges  sometimes  cannot  be  discriminated  within  a  small 
image  sub-area  or,  conversely  that  what  appears  to  be  an  edge 
within  the  sub-area  could  belong  to  a  homogeneously  textured 
domain  of  the  picture.  Increasing  the  size  of  the  sub-area 
apparently  solves  the  problem,  but  is  limited  by  computa¬ 
tional  costs  and  the  complexity  of  large  segment  description. 

It  is  generally  recognized  that  boundary  detection  is 
therefore  best  done  by  the  combination  of  a  relatively  simple 
edge  and/or  line  segment  detector,  followed  by  algorithms 
that  thin  and  link  the  segments  obtained  into  continuous 
boundaries.  Several  fast  numerical  techniques  for  luminance 
edge  extraction  have  been  published,  for  example,  Robert's 
"gradient"  [1 ] ,  Kirsch" s  [2]  ,  Sobel ' s  [3],  Prewitt's  [4], 
Robinson's  [5],  and  the  so-called  "smoothed  gradient"  [6] 
operators.  Comparison  of  the  above  algorithms  reveals  simi¬ 
larities  that  suggest  underlying  general  principles. 

From  these  we  develop  a  set  of  orthogonal  functions 
which  are  closely  related  to  distinctive  image  features.  The 
properties  of  these  functions  suggest  ways  to  minimize  the 
amount  of  computations  as  well  as  an  improved  decision  crite¬ 
rion.  Considerable  improvements  are  obtained  in  terms  of 
boundary  "thickness"  and  sensitivity  to  faint  edges.  The 
edge  detection  methods  are  described  in  Section  3.2. 

The  simplicity  of  edge  data  with  its  low  computer 
memory  requirement  suggests  several  ultra- fast  image  registra 
tion  methods  through  edge  correlation.  Theories  of  edge 


correlation  have  been  investigated  and  several  experiments 
have  been  performed  to  test  their  applicability  to  real 
images.  Section  3.3  describes  edge  registration  as  a  poten 
tial  method  of  image  matching.  Experiments  were  performed 
to  register: 

*  Optical-to-optical  images 

*  Radar-to-radar  images 

*  Radar-to-optical  images 

The  results  of  these  experiments  are  discussed  in 


Section  3.3. 


68 


3.  a..  Fast  Edge  Detection 
3.2.1  Definitions 

The  problem  of  boundary  element  detection  can  be  for- 

2 

mulated  as  follows:  given  a  set  of  n  luminance  samples  from 
an  image  sub-area,  determine  whether  the  sub-area  contains  a 
boundary  element  between  two  regions  of  different  homogeneous 
luminances  (edge) .  It  may  also  be  of  interest  to  determine 
whether  the  area  contains  a  line  or  a  pair  of  degenerate 
edges  enclosing  an  object  too  thin  to  be  resolved.  To  this 
end,  we  define  the  following  mdoels  of  "ideal  boundary  ele¬ 
ments  . 

Consider  an  image  sub-area  A  of  size  n  x  n  sampling 
intervals  shown  in  Figure  31a.  In  the  continuous  image  do¬ 
main,  we  define  an  "ideal  edge  element"  as  a  straight  boun¬ 
dary,  passing  through  the  center  of  A,  and  which  separates 
two  regions  of  different,  constant  luminances  b^  and  b^. 
Adopting  the  convention  b^  >  b^,  the  direction  $  of  the  edge 
element  is  uniquely  determined  with  respect  to  any  arbitrary, 
fixed  direction  as  shown  in  Figure  31b.  The  ideal  edge  ele¬ 
ment  is  characterized  by  its  "magnitude"  =  | ^  -  b^ I  and 
orientation  g,  ^  2-it  . 

Next  we  define  an  "ideal  line  element"  (in  the  continu¬ 
ous  image  domain)  as  a  straight  strip  of  width  approximately 
equal  to  one  sampling  interval,  passing  through  the  center  of 
A,  and  of  different  luminance  b^  than  its  surrounding  b ^  illu¬ 
strated  in  Figure  31c.  The  ideal  line  elements  is  charac¬ 
terized  by  its  "magnitude"  =  ]  b^  -  b 2 1,  its  orientation  4>  ^ , 

Oi.  <f>  and  polarity  sign  (b-L  -  b2)  . 


70 


Finally,  an  "ideal  point"  is  defined  as  a  point  of 
brightness  b-^  different  than  a  constant  brightness  b^  of  the 
surround.  The  ideal  point  is  characterized  by  its  "magni¬ 


tude"  =  | b^  -  b2 J  and  polarity  sign(b^  -  b^)  . 


For  the  discrete  case,  we  define  the  following  notation. 


Consider  the  set  of  n  luminance  samples  b^_.  of  the  image  sub- 

2 


area  as  an  element  of  an  n  -dimensional  vector  space  B.  The 
elements  of  B  can  be  represented  by  a  matrix  B  or  a  column 
vector  b,  for  example  (n=3) 


’bll 

b12 

1 - 

CO 

ja 

B  = 

b21 

b22 

b23 

-b31 

b32 

b33  - 

,  or  b  = 


'1 


(1) 


Finally  we  define  an  inner-  (or  dot)  product  (*,')  on  B  as 


n  n  a 

(B,  C)=  E  S  b..c  or  (b,  c)  =  E  b.c.  . 

•  _  i  •  i  11  LI  "  1  1 


(2) 


3.2.2  Review  of  Previous  Work 

Previous  fast  edge  detection  algorithms  [1  -5]  fall 
essentially  into  two  categories: 

a)  Evaluate  the  maximum  average  gradient  AG  (or 
"smoothed  gradient")  present  in  each  image  sub-area.  The 
average  is  estimated  in  a  direction  perpendicular  to  the 
(unknown)  edge  element  orientation  and  the  maximum  is  approxi¬ 
mately  obtained  by 

AG  «  [(B,  Wx)2+ (B,  W2)2]^  (3) 

where  B  is  the  vector  of  luminance  samples  and  ,  W2  are 
weighting  functions  shown  in  Figure  32  a,  b,  c[l,3,6].  When 
the  average  gradient  exceeds  an  arbitrary  threshold,  the  image 
sub- area  is  considered  to  contain  an  edge  element.  Orienta¬ 
tion  is  then  obtained  approximately  as 

<p=  atan[(B,  Wj)/(B,  W2)].  (4) 

In  order  to  simplify  computations,  the  sum  of  squares 
of  eg.  3  is  sometimes  replaced  by  a  sum  of  absolute  values. 

It  is  oointed  out  that  the  above  measures  are  not  isotropic, 
e.g.  certain  edge  orientations  are  favored  over  other  ones 
[5],  A  pair  of  isotropic  weighting  fuctions  is  shown  in 
Figure  32d. 

b)  The  second  approach  is  to  form  inner  products  of  the 
luminance  vector  B  with  a  set  of  discrete  edge  templates  or 
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masks  T,  of  different  orientations  as  shown  in  Figures  33a,  b, 
c,  d[2],  [4],  [5]),  and  retain  the  largest  value 

max  {  (B,  T\  )  }  . 

When  this  value  exceeds  an  arbitrary  threshold,  the  sub-area 
B  is  considered  to  contain  an  edge  element.  The  direction  is 
approximately  equal  (±tt/4)  to  the  orientation  of  the  template 
giving  the  largest  inner  product. 

This  second  conceot  can  be  immediately  extended  to  line 
and  point  detection  with  the  template  vectors  shown  in 
Figure  34a  and  b. 


I 
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*  3.2.3  Derivation  of  Orthogonal  Feature  Basis 

1  We  now  seek  an  appropriate  basis  for  B.  Because  the 

templates  of  Fig.  33c  and  34a  reoresent  samples  of  ideal  edge 
and  line  elements  positioned  in  eight  equidistant  orienta¬ 
tions,  we  assume  "edge"  and  "line"  subspaces  of  B  spanned  by 
these  vectors,  of  all  possible  orthogonal  bases  for  these 
|  subspaces,  we  choose  the  one  shown  in  Fig.  35  because  of  the 

following  properties:  a)  the  first  pair  of  basis  vectors 
and  VI 2  represents  the  isotropic  smoothed  gradient  weighting 
function.  This  oair,  taken  together  with  the  second  pair 
spans  the  above  "edge"  subspace,  b)  The  second  pair  of  basis 
vectors  and  has  a  distinctive  higher  order  aspect 
(three  zero  crossings  instead  of  one)  and  will  be  shown  to 
contribute  little  to  the  magnitude  of  the  edge  subspace  com¬ 
ponent.  c)  The  "line"  vectors  were  decomposed  into  pair  of 
vectors  W,.,  Wg  with  directional  preference  and  a  pair  ,  Wg 
without  directional  preference.  Note  that  the  point  basis 
vector  of  Fig.  34b  is  equal  to  the  sum  of  the  latter  pair, 
which,  incidentally,  span  all  possible  discrete  realizations 
of  the  discrete  Laplacian  [ 7 1 .  Finally,  the  vector  Wg  was 
added  to  complete  the  basis.  Observe  that  linear  combina¬ 
tions  of  each  pair  of  vectors  produce  similar  distinctive 
patterns,  which  we  call  "average  gradient,"  "ripple,"  "line,” 
and  "Laplacian"  respectively. 

Fig.  36  and  37  illustrate  the  above  discussions.  An 
original  image  of  size  256  x  256  pixels  was  projected  onto 
each  one  of  the  nine  orthogonal  basis  vectors  of  Fig.  35. 
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Since  the  projections  are  bipolar,  a  constant  value  was  added 
for  display,  and  the  images  were  scaled  for  better  display. 

The  complementary  nature  of  (edge,  ripple)  and  (line,  point) 
spaces  is  clearly  visible  and  can  be  attributed  to  the  fact 
that  the  former  basis  vectors  are  odd  with  respect  to  one 
axis  of  symmetry  whereas  the  latter  basis  vectors  are  even 
(see  Figures  33c  and  34a)  . 

In  order  to  reduce  computation,  it  would  be  desirable  to 
reduce  the  dimension  of  the  feature  subspaces.  Figure  38 
shows  the  magnitudes  of  the  projections  onto  "average  gra¬ 
dient,  "  "ripple,"  "line,"  and  "point"  respectively.  While 
the  last  pair  of  projections  appear  similar  (any  line  is  com¬ 
posed  of  points!),  the  "average  gradient"  and  "ripple”  pro¬ 
jections  are  quite  different.  It  is  easy  to  see  that  the 
"ripple"  subspace  adds  little  to  the  "edge"  subspace  and  may 
be  ignored  to  save  computation. 

An  improved  edge  picture  can  be  obtained  by  a  simple  use 
of  the  complementary  property  of  projected  components.  If  we 
subtract  the  edge  magnitude  form  the  sum  of  line  and  point 
magnitudes,  the  resultant  picture  will  have  diminishing  values 
of  those  "edge  points"  with  impure  edgeness,  making  real 
edges  sharp  and  obvious.  Depending  on  the  noise  level  of 
the  original  image  an  alternative  edge  decision  criterion 
could  be  a  relative  edgeness  thresholding  in  stead  of  absolute 
one.  In  this  way  faint  edges  can  be  detected  with  the  same 
opportunities  as  those  of  strong  edges.  However,  care  must 
be  taken  as  deciding  faint  edges  because  they  are  more  sub¬ 
ject  to  noise. 
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Magnitude  in  E-Subspace  Magnitude  in  L-Subspaco 


Magnitude  in  P-Subspace  Magnitude  in  R -Subspace 


Figure  38. 
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Following  this,  the  edge  picture  can  even  be  improved 
by  running  the  convolution  (Fig.  34b)  on  it.  This  'Doint 
operation"  effectively  extracts  peak  values  of  the  source 
array,  making  the  edges  even  thinner  and  prominent.  However, 
this  follow-up  process  will  raise  the  noise  level  somewhat 
and  it  has  to  be  accompanied  by  noise  cleaning.  This  is 
similar  to  peak  detection  in  [9], 


! 
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3.3  Edge  Registration 

3.3.1  Edge  Correlation  Methods 

For  the  purpose  of  fast  feature  mapping  the  edge  infor¬ 
mation  of  both  the  search  region  and  searching  windows  is  ex¬ 
tracted.  Correlation  or  similarity  algorithms  are  then  run 
on  them.  There  are  several  advantages  to  edge  feature  regi¬ 
stration  instead  of  original  scene  magnitude  registration. 
Firstly,  edges  are  more  likely  to  be  sensor  independent 
(optical  vs.  radar).  Secondly,  they  are  invariant  to  illumi¬ 
nation  change  or  even  to  a  local  illumination  reversal. 
Thirdly,  the  storage  of  edge  map  requires  up  to  88%  less 
storage,  because  of  the  low  percentage  of  edge  points  among 
a  picture  in  general.  If  bit  operation  is  available  the 
storage  can  even  go  down  for  several  other  orders.  Finally, 
search  computations  can  be  reduced  to  logical  functions  in¬ 
stead  of  fixed  point  multiplications. 

Let  two  images,  A  the  search  region  and  B  the  window  be 
defined  as  M  x  N  and  m  x  n  pixels  respectively  and  both  are 
represented  by  integers.  The  task  is  to  find  the  sub-area  of 
A  which  registers  or  best  fits  B  (Appendix  A) .  Either  a 
sequential  or  random  search  strategy  can  be  implemented. 
Besides  this,  there  are  several  ways  to  define  similarity  and 
their  computational  requirements  differ  considerably!!  1  0]  . 

The  elements  of  unnormalized  cross  correlation  surface 
R(i,j)  are  defined  as 

n  m 

^  (i/j)  “  ^  (k» £)  •  B  (i+k,  j-htj  (5) 


A  search  for  the  maximum  on  the  surface  is  then  initiated. 


However,  the  maximum  does  not  necessarily  yield  the  registra¬ 
tion  point  or  yields  a  bad  fix.  Therefore,  the  correlation 
must  be  normalized  as 

n  m 


£  £  A  (k,  i) .  B  (i+k,  )  +  D 


n  m 

'  n  m 

\-l 

B2(i+k,;+t); 

a 

(6) 


The  resultant  fix  is  much  better  than  the  unnormalized  cor¬ 
relation  and  the  normalized  correlation  does  give  excellent 
results  as  to  the  self  correlation  of  both  optical  and  radar 
scenes.  It  is  sensitive  to  background  level  and  it  requires 
extensive  computations. 

As  a  remedy,  the  followino  new  similarity  measures  are 
proposed.  Further  details  of  which  are  described  in  Appendix 
B.  These  measures  represent  the  edge  pictures  with  binary 
numbers  0's  and  l's  for  edge  points  and  non-edge  points  re¬ 
spectively.  Define  a  logical  similarity  measure  as 


n  m 

£  £  A  (k,-C)  -  B  fk+i.-C+i) 


k=n=i 

■  n  m 

n  m 

h 

£  ^  A  (k,t)  X 

£  £  B  (k+i,  £+j) 

.k  =  t  l  =  \ 

k=l  t=l 

(7) 


formula  is  for  use  in  the  case  where  only  edge  points  are 
to  be  considered  and  can  be  further  approximated  as 
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Alternatively,  if  all  points  are  to  be  considered  (edge 
and  non-edge  points)  we  define  the  NOT-EX-OR  correlation 
measure  as 


EXOR 


— 1 £  £  _  - 

m  *  n  ,k=i  1  =  1  A  (k’l>  e*-or  b  fk+i, /„+,*) 


Note  that  (i,  j),  RAND(i,j)  and  Rjjxqr  have  values  between 
0  and  1;  they  are  normalized  quantitives.  A  comparison  of 
the  different  correlation  computations  is  shown  in  Figure  39. 
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3.3.2  Experimental  Results 

A  sequence  of  experiments  were  performed  using  the 
edge  correlator  method  (Eq.  (8))  to  register: 

•  Optical-to-optical  images 

•  Radar- to-radar  images 

•  Radar-to-ODtical  images 

Figure  40  shows  the  edge  pictures  extracted  from  the 
optical  and  radar  images.  The  size  of  both  images  is  180  x 
180.  A  window  of  size  64  x  64  located  at  the  coordinates  of 
(57,  57)  was  selected  from  the  optical  image  as  shown  by 
solid  lines  in  Figure  40(a).  This  window  was  correlated  with 
the  180  x  180  optical  image  at  every  fourth  location  to  pro¬ 
duce  a  correlation  map  of  size  30  x  30.  The  result  is  shown 
in  Figure  41a.  Similarly,  the  radar- to-radar  correlation 
was  performed  and  the  results  are  in  Figure  41b.  Both 
figures  show  sharp  correlation  peaks  located  at  the  points 
of  true  registrations,  thus  providing  a  high  confidence  level 
that  the  correct  matches  have  been  found.  These  results 
indicated  that  when  two  images  are  relatively  free  from  geo¬ 
metric  and  sensor  distortions,  the  implementation  of  edge 
correlator  using  the  fast  algorithm  with  the  AND  operator 
as  in  equation  (8)  is  especially  attractive  in  real  time 
computation.  It  reduces  computations  and  storage  require¬ 
ments  since  the  fixed  point  multiplications  have  been  replaced 
by  the  simple  AND  operations. 

To  see  how  the  edge  correlator  can  be  applied  to  radar- 
to-optical  registration,  two  windows  of  different  sizes  and 
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Figure  40.  Optical  and  radar  edge  pictures  1180  x  180) 


a.  Optical-to-optical  correlation 


b.  Radar-to-radar  correlation 


Figure  41. 
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locations  were  chosen  from  optical  images.  The  two  windows 
were  then  correlated  with  the  radar  image. 

The  window  sizes  and  locations  as  shown  in  Figures  40 
and  43  can  be  summarized  as  follows 

Window  No.  Window  Size  Location  Reference  Size 

1  (64  x  64)  (57,  57)  (180  x  180) 

2  (32  x  64)  (49,  21)  (  90  x  122) 

Window  No.  1  was  correlated  with  the  reference  image  at 

every  4th  location.  Window  No.  2  was  correlated  with  the 
reference  image  at  every  2nd  location  Figure  42  shows  the 
correlation  of  window  No.  1.  For  clarity  a  threshold  was 
selected  so  that  only  those  points  with  relatively  high  cor¬ 
relation  were  shown.  Similarly  the  correlations  of  window 
No.  2  are  shown  in  Figure  44.  In  the  case  of  window  No.  1, 
althouah  a  relatively  high  Deak  was  found  to  be  located  at 
the  true  registration  location,  other  peaks  of  equal  or 
higher  amplitude  were  also  obtained.  This  provided  only  a 
relatively  low  confidence  in  locatim  the  true  registration 
point.  In  the  case  of  window  To.  2,  a  high  peak  was  obtained 
at  a  location  which  is  four  pixels  away  from  the  true  re¬ 
gistration  location.  An  examination  of  the  optical  window 
and  the  reference  radar  map  indicates  that  a  slight  geometric 
misregistration  exists  between  the  two  images.  Thus,  the 
correct  matching  location  was  located  more  accurately  by  the 
automated  technique. 
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3.4  Conclusions 

Several  significant  techniques  on  edge  extraction  and 
registration  have  been  developed  in  this  section.  These  are: 

1.  A  set  of  edge  operators  extracted  from  the  derived 
orthogonal  basis  were  identified.  These  operators  encom¬ 
pass  nearly  all  of  the  previously  derived  operators  by  other 
researchers  on  a  3  by  3  image  window.  Therefore,  the  new 
edge  operators  can  be  applied  to  extract  edges  as  well  as 
other  important  image  features  such  as  points  and  lines. 

2.  A  correlator  with  AND-operations  was  developed.  This 
correlator  is  far  more  efficient  than  the  existing  correla¬ 
tors  in  terms  of  computational  efficiencies  and  storage  re¬ 
quirements.  Theoretical  and  experimental  results  derived 
from  this  section  indicate  that  edge  feature  can  be  used  as 

a  useful  tool  in  image  registration.  These  results  indicate 
that  the  edge  correlator  is  efficient  and  works  very  well  on 
most  regions  of  a  scene  in  which  only  a  few  salient  edges  are 
present  and  no  significant  geometric  or  sensor  distortion  is 
present.  Whenever,  geometric  or  sensor  distortions  are  pre¬ 
sent,  the  edge  correlator  performance  is  heavily  dependent 
upon  the  scene  content  as  well  as  the  degree  of  distortion. 
Further  research  is  needed  to  predict  the  degree  of  sensiti¬ 
vity  to  these  distortions  as  well  as  on  other  edge  detection 
operators  and  techniques. 
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CONVENTIONAL  SIMILARITY  MEASURES 
FOR  CORRELATION  SEARCH 


-  Search  region  A  is  M»N  pixels 

-  Window  scene  B  is  m-  n  pixels 

-  A  and  B  are  represented  by  integers 


Cross  -  Correlation 


n  m 


R  (1>j)  ■  EE 


k=l  -6=1 


A  (k,  l)  •  B  (i  +  k,  j  +*,) 


Problems:  R(i,  j)  is  a  function  of  energy  content  of  B  within 
window  at  position  (i,  j).  Fix  is  very  bad  (see  2-D 
plot) 


Normalized  Cross-Correlation 


RN(i,j)= 


A  (k,  l)  B  (i  +k,  j  +  -t ) 


231^  [A(k,t)]2  •  2353  [  3(i+k,j+t)]2] 
L 


* 


Remarks:  Fix  is  greatly  improved  (see  2-D  plot) 
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Appendix  B 

THE  LOGICAL  NORMALIZED  CROSS- CORRELATION  MEASURE 

A  NEW  SIMILARITY  MEASURE  FOR 
CORRELATION  SEARCH  OF  EDGE  MAPS 


Problem:  The  conventional  cross-correlation  or  normalized  cross¬ 

correlation  measures  cannot  be  used  for  binary  representation; 
bipolar  representation  must  be  used  to  account  for  points  where  A=0, 
B=1  or  A=l,  B=0.  Failure  to  do  so  would  bring  in  only  points  where 


A=B=1. 


Binary  representation  of  edges 


Bipolar  representation  of  edges 


1  n  n  n 

0 — 1  -f— J/bUl — r* 

edge  points  ?c>an 

line 

w-.n-_-_-ji-.-n-— 1 L 


This  is  achieved  at  great  computational  savings  by  the  LOGICAL 
NORMALIZED  CROSS- CORRELATION  MEASURE 


-  The  sums  are  simple  incremental  counters. 

-  The  integer  multiplications  of  the  numerator  are  replaced  by  a 

1  bit  EX  -OR  function. 

-  The  integer  squares  of  the  denominator  are  eliminated. 

-  Sequential  similarity  search  can  be  implemented  by  randomizing 

the  memory  addresses  corresponding  to  the  indices  k,  l, 
k+i,  -t+j. 
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SECTION  IV. 


INVARIANT  MEASUREMENTS 


4.1  INTRODUCTION 

The  theory  of  algebraic  invariants  from  which  the  various 
facets  of  invariants  developed  is  based  upon  a  simple  idea 
which  first  appeared  in  Lagrange's  papers  and  was  subsequen¬ 
tly  developed  by  Gauss  [1,2]. 

This  simple  idea  was:  given  a  quadratic  equation 
2 

ax  +  2bx  +  c  =  0 
lets  us  transform  x  to  y  by: 


=  p*+q 


or 


rx+s 

then  transformed  equation  will  be 


q-sy 

ry-p 


2 

Ay  +  2By  +  C  =  0  where 
A  =  as2  +  cr2  -  2bsr 

B  =  -aqs  +  b(qr+sp)  -  cpr 

2  2 
C  =  aq  -  2bqp  -  cp 

then  let  us  form  the  discriminant: 

B 2  -  AC  =  (ps-qr)2(b2  -  ac) . 

Hence  the  discriminant  of  the  equation  in  y  is  equal  to  the 

2 

discriminant  of  the  x  equation,  times  the  factor  (ps-qr) 

which  depends  only  upon  the  coefficients  p,  q,  r,  s  in  the 

transformations  y  =  . 

rx+s 

Boole  was  the  first  to  observe  that  for  every  equation 
the  discriminant  remains  unchanged  (except  for  a  factor)  if 
x  is  transformed  to  y  by  some  transformations. 


i 


A  uniform  method  which  would  give  all  the  invariant 
expressions  was  given  by  Cayley  in  1945  by  his  memoir,  "On 
the  Theory  of  Linear  Transformations."  Then  Sylvester  joined 
him  in  the  construction  of  the  general  theory  of  algebraic 
invariants . 

The  development  of  invariant  algebra  continued  through 
the  first  half  of  the  present  century  by  efforts  of  E.  Elliot, 
Salmon,  Grace  and  Young,  et.  al. 

Algebraic  invariant  theory  contributes  the  derivation 
of  algebraic  expressions  which  will  remain  invariant  under 
certain  transformations.  The  first  use  of  algebraic  invari¬ 
ants  and  moment  invariant  measurements  for  the  purpose  of 
picture  recognition  was  advanced  bv  Hu  [2] .  He  used  spatial 
moments  as  a  means  of  representing  characters  and  normalized 
them  by  developing  a  theorem  relating  momement  measurements 
to  algebraic  invariants.  Hall,  et  al .  [3,4]  has  used  moment 
invariants  to  describe  both  the  spatial  distributions  of 
objects  and  edge  structures  of  texture  patterns  for  chest  X- 
ray  recognition. 

A  major  problem  encountered  in  the  computer  representa¬ 
tion  and  measurement  of  a  scene  is  the  fact  that  by  changing 
the  viewing  angle,  scale,  or  rotation,  most  numerical  measure¬ 
ments  will  change  substantially.  Therefore,  it  is  desirable 
to  consider  normalization  of  the  derived  measurements  in 
such  a  way  that  invariant  properties  be  established. 


Any  discrete  two  dimensional  scene  may  be  represented 
by  f(x^,y^);  i  =  1,2,...,N,  j  =  Certain  measure¬ 

ments  from  this  two  dimensional  array  may  be  normalized  by 
the  utilization  of  algebraic  invariant  techniques. 

The  purpose  of  this  section  is  to  explore  more  deeply 
the  use  of  moment  invariants  for  the  analysis  and  represen¬ 
tation  of  complex  scenes,  present  an  analysis  of  a  variance 
of  the  "invariant  measurements"  when  made  on  computer  images, 
and  present  pictorial  examples  for  scene  analysis. 
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4.1.2  QUANTICS  AND  INVARIANTS 
Quantic 

A  function  of  several  variables  x, y, z, . . .which  is  ra¬ 
tional,  integral  and  homogeneous  is  called  a  quantic. 

If  there  are  only  two  variables,  such  a  function  is 
called  a  binary  guantic,  and  for  a  qeneral  g  variable  we 
call  it  q-ary  quantic. 

The  order  of  a  quantic  is  the  degree  of  x,y, .... 

The  quantics  of  1st,  2nd,  34d,...,  pth  are  called 
linear,  quandratic,  cubic,...,  p-ics  respectively.  The  binary 
p-ic  will  be  denoted  by  the  form:  a^x^*  +  • • 

arx*',-r+'*'yr  which  is  usually  shown  in  short  as: 

p 

(aQ,a1,a2, . . .,ap)  (x,y)  . 

Note:  When  we  talk  about  the  coefficient  of  a  binary  quantic 

we  mean  a.,d.  ,  ..  *,a  and  not  a .  ,  pa la  ,...  . 

01  p  0^0'  \  r  I  r 


p-r+1 

r 


102 


Invariants 

An  invariant  of  a  single  quantic  is  a  function  of 
the  coefficients  of  that  quantic  which/  after  a  linear  trans¬ 
formation,  remains  constant  except  for  a  multiplicative 
factor  which  is  a  function  only  of  the  coefficients  of  the 
transformation . 

Ex.  Consider  a  binary  quantic  of  order  p: 

Q  =  (ao,ai'**-»ap)(x»  y)p 

in  which  we  transform  x,y  by 


then  the  same  quantic  in  the  transformed  domain  would  be 


Q  = 


a^Kx1, 


where  the  values  of  a^,a^,...,a^  are  evaluated  by  the  iden¬ 
tity: 


(aQ*  aj*  •  •  • »  ap)(x»  y)P  =  (ag,  a  ,  a  )(x  ,  y  )P  . 


Now  we  say  f(ag,...,a  )  is  invariant  under  the  above  trans¬ 
formation  if  the  following  identity  holds: 


V  =  ^^11'  ^12’  ^21'  ^22^aQ'  ’ 


•V 
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where  0  is  a  power  of  the  determinant  of  the  transformation: 

A  “ 

p  =  A 

A  being  the  determinant  of  the  transformation,  i.e. 

A  =  l ^5. 22" a ^2^ 21  ^or  t*le  case  “  =  0  we  bave  absolute 
invariants . 

Two  sets  of  (x,y, . . .)  and  (x  1 ,y' , . . .)  are  called  co¬ 
gradient  if  a  linear  transformation  on  one  set  compells  the 
other  set  to  the  same  linear  transformation. 

Two  sets  (x^  >*2'  •  •  •  ,xo>  an<^  >  ¥2  ’  *  ‘  ’ '  are  ca-*-le<^ 
contragradient  if  we  transform  the  first  the  first  set  by 

x’  =  Lx 

where 


L  rjJ 


det  X,  4  0 


r-  x; 


x'  = 


UX  j 


X  = 


r  X ,  1 


then  the  transformation  for  the  second  set  has  to  be 

Z*  =  h~Ty 
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where 


-T  . 

L  is  the  transpose  of  the  inverse  of  L 


Theorem:  In  contragradient  quantities  x  y  remains  invariant. 


Proof:  x'  =  Lx 


i'  =  ifTx 

,  T  ,  T_  T  -T  T 

x  y  =  x  L  L  y  =  x  y 


4.1.4.  Invariants  of  Orthogonal  Transformations 
Suppose  we  have  a  binary  quantic  of  order  p: 

(ap0'ap-i,i'---'a0p<!!'y)P 

consider  the  transformation 


cos  0  -sin  0\  /  x' 
sin  0  cos  0  /  \  y 


Now  we  transform  (x,y)  and  (x',y')  as  follows 


then  we  would  have 


Then  the  invariants  of  which  there  are  p+1  linearly  indepen¬ 
dent  ones,  are  obtained  from  the  following  equations. 


I'  =  eLP0T 
pO  e  pO 


r'  =  ei(p-2)0 

P-1.1  P-1.1 


r  /  =  -i(p-2)GT 

h,p-l  6  h,p-l 

I  /  =  --1P0t 

Op  *0p  * 


From  the  first  two  identities  we  have 


Vo  *pO‘l(l)mp-».r  (2  )ap-2f  z 


rp-l,  l  =  (ap0+ap-2,  Z]  '  l(P'2)(ap-l,  l  +  ap-3,  3} 
+  ...  +  (-i)P"2(a2>p.2+a0p) 


^p-2.,  2  =  (ap0+2ap-2,  2+ap-4,4)  ‘  l(p_4){ap-l,  l+2ap-3,  3+ap-5,  5* 
+  ...  +(-i)P“4(a4>p.4+2a2p2+a0p) 


I  =  [(a  a  ;a  0  0  )(l,  1) r; 

p-r,  r  1  pO  p-2,2  p-2r,  2r"  '  ’  * 

(ap-l,  lJap-3,  3’  *  *  *  ;ap-2r-l,2r+l)(1'  1)r''*** 


(a2r,p-2r;a2r+2,p-2r-2;"'  ;a0p,(1'  -i)P"2rp-2r>  0 
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4 .2  Perceptual  Invariants 

In  a  discussion  of  shape  recognition,  Deutsch[3]  de¬ 
scribed  four  invariance  properties  of  human  perception. 

1.  Shapes  may  be  recognized  independent  of  their 
location  in  the  visual  field.  It  is  not  necessary  to 
fixate  on  the  center  of  a  figure  in  order  to  recog¬ 
nize  it  nor  need  the  eyes  be  moved  around  the  con¬ 
tours  of  a  figure.  Thus,  shape  recognition  is  in¬ 
variant  to  translation  provided  the  object  remains 
in  the  visual  field. 


2.  Recognition  can  be  effected  independently  of 
the  angle  of  inclination  of  a  figure  in  the  visual 
field.  The  angle  refers  to  the  angular  orientation 
of  the  figure  in  two  dimensions  although  some  in¬ 
variance  to  the  depth  angle  is  also  observed. 
Therefore,  the  recognition  of  a  figure  is  invariant 
to  rotation  in  two  dimensions,  although  the  angle  of 
rotation  may  also  be  recognized. 

3.  The  size  of  a  figure  does  not  interfere  with 
the  recognition  of  its  shape  provided  the  entire 
figure  is  within  the  visual  field.  Thus,  visual 
recognition  exhibits  an  invariance  to  scale. 

4.  Mirror  images  appear  alike.  Both  humans  and 
animals  tend  to  confuse  these.  This  observation 
rules  out  a  template  matching  theory  of  shape  re¬ 
cognition  because  a  template  superposition  cannot 
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take  place  for  mirror  images  and  thus  the  confu¬ 
sion  cannot  be  predicted. 


Deutsch  also  points  out  that  certain  primitive  organisms 
find  it  more  difficult  to  distinguish  between  squares  and 
circles  than  between  rectangles  and  squares, which  casts 
some  doubt  on  recognition  theories  based  on  angular  properties 
of  figures.  He  also  argues  that  the  invariance  properties 
are  common  to  all  parts  of  the  visual  cortex,  which  rules  out 
a  scanning  process  for  shape  recognition.  To  illustrate  that 
the  invariance  properties  could  be  realized  easily,  he  pro¬ 
posed  the  following  system. 

The  receptor  surface  is  connected  to  a  two  dimensional 
coding  array  which  is  excited  by  a  contour  falling  on  the 
corresponding  receptors.  For  example,  if  a  spherical  object 
is  focused  on  the  retina,  the  circular  contour  is  projected 
on  the  coding  array.  The  processing  of  the  edge  contour  is 
described  by  the  following  four  rules. 

(i)  Every  figure  excites  the  retina  only  by  its 
contour  i.e.,  only  the  contour  of  the  figure 
will  appear  on  the  to  w  dimensional  array. 

(ii)  Every  array  point  will  transmit  a  pulse  down 
the  common  cable  as  soon  as  an  excitation 
reaches  it  from  the  retina.  The  integration 
of  tlV'Se  pulses  gives  a  measure  of  the  total 
number  of  edge  points. 

(iii)  Each  contour  point  on  the  array  will  also  ex¬ 
cite  all  the  units  which  are  at  right  angles 
to  it  on  the  array.  These  excitations  will 
not  pass  to  the  common  cable  but  only  to  the 
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neighboring  cells.  It  is  also  assumed  that 
these  lateral  excitations  move  at  the  same 
speed  in  all  directions. 

(iv)  As  such  lateral  excitation  from  a  point  in  a 
contour  advances  another  message  will  be 
sent  down  the  common  cable  as  soon  as  it  co¬ 
incides  with  another  contour  point. 

The  recognition  of  shape  is  made  from  the  signals  on 
the  common  cable.  For  each  pattern,  two  sets  of  pulses, 
separate  in  time,  are  transmitted.  The  first  set  relates  to 
the  number  of  contour  points.  The  second  provides  shape  in¬ 
formation.  This  system  can  be  shown  to  be  invariant  to 
translation,  rotation,  scale  and  mirror  image  as  previously 
described.  The  invariance  to  translation,  rotation  and  mirror 
image  transformations  would  affect  neither  the  number  of  con¬ 
tour  points  nor  the  distance  between  edges  in  the  figure. 

The  size  invariance  may  be  obtained  by  normalizing  the  se¬ 
cond  signal  by  the  first.  That  is,  the  lateral  excitation 
signal  would  be  used  to  normalize  the  occurrence  times  of  the 
second  signal . 

It  is  interesting  to  note  that  the  proposed  system 
would  be  more  easily  confused  in  differentiating  circles  and 
squares  than  rectangles  and  squares. 

Both  circles  and  squares  would  produce  two  sharp 
pulses  which  could  be  of  equal  magnitude  if  the  same  number 
of  edge  points  were  detected.  However,  the  rectangle  would 
j  produce  three  pulses  and  could  easily  be  distinguished  from 

1 

i 

i 

i 


112 


either  a  square  or  circle. 

The  limitations  of  this  simple  model  are  many;  however, 
it  does  illustrate  a  logical  approach  to  recognition.  First, 
the  desired  recognition  properties  are  determined,  then  a 
system  which  satisfies  these  properties  is  designed.  The 
particular  system  described  was  devoted  to  shape  recognition. 
Other  systems  need  to  be  developed  which  are  invariant  to 
intensity  changes  and  more  complex  images.  The  power  and 
simplicity  of  Deutsch's  approach  was  a  major  contribution 
from  psychology  to  the  understanding  of  recognition  systems. 
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4.3.1  Moment  Invariants 

Any  two-dimensional  pattern  can  be  represented  by  a 
picture  function,  f(x,y),  with  respect  to  a  pair  of  axes 
fixed  in  the  visual  field.  It  is  also  true  that  by  means  of 
a  uniqueness  theorem[3],  two  dimensional  moments,  m  , 

fr  Si 

which  are  defined  as: 


m 


pq 


x^y^£(x,  y)dxdy 


p,  q  =  0,  1,2,.. 


can  be  used  to  uniquely  represent  f (x,y) .  Hence  one  may  use 
m  as  a  means  of  representing  any  two-dimensional  pattern. 

The  computations  of  m  consist  of  multiplying  the 
functions  f(x,y)  by  a  monomial  x^y^  and  integrating  the  re¬ 
sult.  The  monomials  of  order  3  or  less,  x^y^,  x^y^,x^y^, 

1011. 2  02130  .12  .  .c 

xy  ,xy  ,xy  ,xy  ,xy  ,  and  x  y  are  shown  m  Figure  45. 

The  moments  of  order  p+q  may  also  be  interpreted  as  the 
response  of  an  imaging  system  with  the  transfer  function  of, 
x^V5,  and  the  input,  f(x,y).  Low  order  moments  have  intui¬ 
tive  relations  to  objects.  For  example,  m^g  is  related  to 
mass  and  m^^  to  center  of  mass  and  m.^,  mQ2'  an(^  m20  to 

principal  directions. 

Another  useful  concept  is  the  moment  generating  func¬ 
tions  which  is  defined  as: 


r"  r°* 

M(u,  v)  =  J  J  exp[^x  +  vy]f(x,  y)dxdy  . 
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If  u  and  v  are  considered  as  complex  variables,  this  expres¬ 
sion  is  a  two  sided  Laplace  transform.  For  the  invariant 
development  both  p  and  v  are  assumed  to  be  real.  This  func¬ 
tion  can  also  be  written  in  the  form  of 


,  -  Z  (u)P  (V)q 

M,U’Vl  =  p5o  q?0  PO  Pi  P! 


where  the  exponential  has  been  expanded  in  its  Taylor  series 
equivalent,  assuming  that  moments  of  all  orders  exist.  This 
equation  shows  that  the  moments  may  be  determined  from  the 
derivatives  of  the  moment  generating  functions  evaluated  at 
the  origin. 

Central  mements  are  defined  as 


«>  • 

Upq=j  J  _(x-£)P(y-y)qp(x,  y)dxdy 


where 


Central  moments  may  be  easily  shown  to  be  invariant  under 
translation  and  from  here  for  the  sake  of  simplicity  we  will 
assume  to  be  using  central  moments. 

To  relate  the  moments  to  algebraic  invariants,  one  may 
first  expand  the  exponential  term  in  the  moment  generating 
'unction  to  obtain 


I 
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M(u,  v)  =  f  r  £  —  (ox  +  vy)Pp(x,  y)dxdy. 

J_oo  p  =  0 

Now  after  using  the  binomial  expansion  and  carrying  out  the 
integration: 


Mi-'=p^v,C)(p;r)(> 


The  functions  of  the  form 


r?oWOC)(?; 


are  called  qualtics  in  the  study  of  algebraic  invariants[l ] 
and  play  a  major  role  in  the  computation  of  moment  invariants. 
The  fundamental  definitions  of  this  theory  are  described  in 
Section  4.2.1. 

4.3.2  Fundamental  Theorem  of  Moment'-  Invariants 
If  the  binary  p-ic  has  an  invariant 

£<a  pO’  *  *  *  *  a0p*  =  ^ap0'  *  *  * '  a0p* 

then  the  moment  of  order  p  has  an  algebraic  invariant 


WQp^  ~  !JU  *  *  * 


»0p> 
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where  J  is  the  Jacobian  of  the  transformation. 

The  importance  of  this  theorem  lies  in  the  fact  one  may 
find  an  invariant  function  of  moments  once  we  have  a  corre¬ 
sponding  algebraic  invariant  function. 

A  point  which  should  be  emphasized  is  the  generality  in¬ 
volved  in  linear  transformations.  The  only  restriction  was 
A  ^  0.  Hence  we  can  use  rotation,  reflection,  magnitude 
change,  and  correspondingly  make  our  normalization  with  re¬ 
spect  to  magnitude,  rotation,  orientation,  etc. 

Example:  Similitude  Moment  Invariants 

Consider  the  transformation 


a  constant  . 


we  can  write  an  algebraic  invariant  simply  as 


a 


P+qa 


pq 


where  a  B  is  the  determinant, 
the  moment  invariant  is: 


By  using  the  fundamental  theorem 


p+q+2p 


pq  ' 


(l) 


Since  J  =  a  . 

For  zero  order  (p+1  =  0)  we  have: 


u’  =  a2y  (2) 

now  by  eliminating  a  between  (2)  and  (1)  be  have  the  following 
absolute  moment  invariants: 


f 

I 


I 

! 


! 

I 


1 


I 


t 


I 
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4.3.3  Derivation  of  Moment  Invariants 

In  the  following  we  attempt  to  use  the  concept  and  tools 
of  invariant  algebra  which  is  discussed  in  Section  4.1  to 
derive  the  seven  invariant  moments  which  are  simultaneously 
invariant  with  respect  to  size,  orientation  and  location  of 
the  object  in  the  scene.  The  detail  theory  behind  these  deri¬ 
vations  is  given  in  Section  4.1.  So  in  the  rest  of  this 
sub  sect'.on  we  proceed  to  use  those  tools  in  calculating 
the  invariant  moments . 


Wl0  =  2S(x-i)1(y-y)0f(x,  y) 


10 


x  y 


=  m 


=  0 


m  i  o 


^11  "  SSiX-xlV-y)1^,  y) 

x  y 


=  m 


11 


miomoi 


m 


00 


^20*  SS(x-x)2(y-y)°f(X/  y) 
x  y 


=  m 


20  m 


2mio2.  +  miQ2 


00 


m 


=  m 


m!02 


00 


20  m 


00 


U02=SS(x-x)  (y-y)  f(x,  y) 
x  y 


m 


moi2 


02  m 


00 
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U30  =  ££<x-x)3(y-y)°f(x,  y) 
x  y 

=  m3()  -  3£m20  +  2m1()i2 

wi2  =  2^(x-x)(y-y)2f(x,  y) 

x  y 

=  m12  -  2ymj  ^  -  xmQ2  +  2y  m1Q 

^21  =  SS(x-x)2(y-y)1f(x,  y) 
x  y 

=  m21  -  2“xmn-  ^n2Q  +  2x2m01 

|i03  =  2Z)(x-x)(y-y)3f(x,  y) 
x  y 

=  m03-  3ym02+  2ym01. 

In  summation: 


^0  =  m00 ' 

UH=mH-ymio 

•l 

o 

ii 

o 

H 

U30=  m30"  3*m20+  2m10i2 

Uoi  =  °» 

U12  =  m12"  2^mll"xm02+2^2ml0 

u20  =  m20  "  *ml0 

w21=m2r  2xmir  ym20+2x2m01 

»02=  moz-ymoi 

wo3 =  mo3"  3ymo2+  2y2moi  • 

The  normalized  central  moments, 


denoted  by 


pq 


are 


defined  as: 
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M 


£SL 


pq  (*?*>) 


u 


00 


p+q  =  2,  3, . . . 


The  normalized  central  moments  are  invariant  to  size 
changes  as  well  as  translation. 

From  the  second  and  third  order  moments  a  last  set  of 
measurements  can  be  derived  which  are  invariant  to  proper 
and  improper  orthogonal  rotation,  and  mirror  images  as  well 
as  invariant  to  translation,  size  change  . 

cpl=  T120+  ^02 

^2=  (t^2o"  T>o2)2+4Tin 
cp3=  (n30-3o12>2+(3o21+Ti03)2 

<P4=  (ri30+r)12)Z+(ri21+Ti03)Z 

2  2 

'Ps  =  ^30“  3r'l2^T130+T112^^r,30+r>12^  “  21+  ^03^  ^ 

+  (3^2r^03){^21  +  T103)r3(Tl30+T1l2)2-(Tl21+^03)21 

2  2 
^r520”'T102^^T130+T112^  ”  ^ 21 +  ^ 

+  4t^  Xl(Tl30+r»12)(Tl21+'n03)  * 

A  seventh  "metric"  can  be  added  which  will  change  sign 
under  improper  rotation. 


4.3.4  Numerical  Computation  of  Moment  Invariants 

In  the  previous  section,  the  moments  were  computed  for 
a  continuous  function  of  continuous  variables.  For  com¬ 
puter  implementation,  the  moments  must  be  computed  for  a 
discrete  function  of  discrete  variables.  Also,  the  range 
of  the  moments  must  be  considered.  In  this  section,  the 
computer  representation  and  numerical  computation  of  moment 
invariants  will  be  described. 

Any  finite  set  of  moments  may  be  computed  and  stored  as 
a  matrix.  For  example,  the  matrix,  M,  of  non-central  moments 
with  p,  q  less  than  or  equal  to  3  is  given  by: 


. 

3 

o 

o 

3 

O 

m02 

m03 

M  = 

m10 

mn 

™12 

ml3 

m20x 

s 

'm21 

m22 

m23 

<;'m30 

m31 

m32 

m33 

If  an  order,  n  =  p+1,  is  selected,  and  all  moments  with  n 
less  than  or  equal  to  fixed  value  is  selected,  then  an  upper 
triangle  matrix  is  produced.  This  condition  is  illustrated 


T 


\ 

\ 
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by  the  elements  above  the  line  in  the  matrix,  M.  Similar 
matrices  may  be  defined  for  central  and  normalized  moments. 

An  important  question  for  the  computer  implementation 
of  moment  invariants  is  the  range  of  the  moments  since  this 
determines  the  computer  storage  requirements.  The  maximum 
value  of  the  noncentral  moments  is  easily  determined  by  the 
following  procedure.  Let  the  objective  function,  J(x,y),  be 
given  by 

J(x,y)  =  mpq. 

To  maximize  J  subject  to  the  constraint 

If f (x, y) dxdy  =  1 

one  may  use  the  method  of  Lagrangian  multipliers: 

H  =  J  +  X  (jjf  (x,y)dxdy  -  1). 

Equating  the  partial  derivatives  of  H  with  respect  to  x,y  and 
X  to  zero,  give  the  conditions 

pf (x, y)  +  xfx (x, y )  =  0 

qf (x , y)  +  yfy(x,y)  =  0. 

Thus,  the  maximum  value  of  mpq  may  be  expressed  as 
Jmax  =  //  ^  "  'r 

X  y 

This  expression  may  be  used  to  determine  the  maximum  value  of 
a  known  function,  f(x,y).  For  example,  suppose  that 

f(x,y)  =  K  exp{-x-y}. 
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Application  of  the  above  formula  gives 


JMAX  =  PP<?q- 


Although  the  previous  expression  gives  the  maximum  value 
of  the  moments,  it  depends  on  the  function.  It  is  sufficient 
for  computer  implementation  to  use  upper  bounds  for  the 
maximum  value.  For  a  finite  rectangular  region,  easily 
computed  bounds  may  be  based  on  the  following  relationship: 


m 


pq 

which  gives 


a  c 


m 


=  J  J  xPyqp(x,y)dxdy  sfnaxJ  J  xPyqdxdy 

3.  C 


s  f _ fbp+1-ap+1  -j^dq+1-cq+1 


q+1 


]• 


pq  max  L  p+1 

A  similar  upper  bound  for  the  central  moments  is  easily  shown 
to  be 


u  */  rtb-a^^IRd.c)^1-] 

Kpq  “  max  L  p+1  -1 q+1  J  * 


Lower  bounds  may  be  computed  in  a  similar  manner. 

The  computed  bounds  on  the  moments  may  be  used  to  scale 
the  computer  operations,  select  register  lengths,  and  norma¬ 
lize  computed  feature  values. 

Another  factor  to  be  considered  in  the  numerical  compu¬ 
tation  of  moment  invariants  is  that  since  the  invariance 
properties  were  proven  only  for  continuous  functions,  certain 
variations  should  be  expected  in  the  discrete  case.  The 
following  development  shows  that  these  variations  do  occur 
but  can  be  controlled  by  careful  techniques. 


Lemma  1.  The  moments  are  invariant  under  translation  of 
coordinates  provided  xQ  =  kAx,  yQ  =  £ Ay. 

Proof:  Suppose  the  density  with  respect  to  coordinates  x,y  is 

f(x,y).  Let  x,*y'  be  a  translation  of  x,y: 


X'  =  X  -  X, 


y'  =  y  -  y„ 


for  some  x  ,  y  .  Then  with  respect  to  the  x',  y'  coordinates, 
o  o 

the  density  must  be  f'(x',y')  =  f(x'+xo,  y'+yo).  With  re¬ 
spect  to  the  x,y  coordinates,  the  p-q*"*1  central  moment  is 

u*.y  =  EE  lrm-?y(* „-rm)&x  iy 

(Vym)Er 


with  respect  to  the  x',y'  coordinates,  the  D-q  central 
moment  is 

Mx.  y.lf'p'’>=EE  (x'„'x'  iVm  f  P  (x'n'  ym)AX' 4y' 

’  mu'  '  '  ' 

(X^y'm)er 


where  of  course.  Ax’  =  Ax,  Ay’  =  Ay i 

Note  that  we  can  index  x'  such  that  x '  =  x  -  x 

n  n  n  o 

iff  xq  =  kAx  for  some  inteqer  K.  Similarly  for  y^.  Hence, 
provided  that  xQ  =  kAx,  yQ  =  £Ay,  then  it  is  readily  verified 
that 


X  '  =  x  -  X 


y '  =  y  -  y,. 


XA  -  X'  =  xn  -  x  ,  y;  -  y  =  ym  -  y 


mu*-,** wm 


n*.  v.  tf'.P.q)  =  uv  v  (f.p.q)  s 

A  > /  y 


i.e.,  provided  that  xq  =  kAx,  y^  =  £  Ay,  then  ( f,p,q ) 

is  invariant  undre  translation  of  coordinates. 

Note:  xq  =  kAx,  yQ  =  £  Ay  is  required  if  (1) ,  (2) ,  (3) , 
(4)  are  to  hold  exactly.  If  xq  /  kAx  and/or  yQ  /  £  Ay, 
then  (1),  (2),  (3),  (4)  become  approximation  and  the  con¬ 

clusion  is  only  aoproximately  valid. 

Lemma  2.  The  moments  are  aoproximately  invariant  under  simi¬ 
larity  transformation. 

Proof:  We  write  the  moments  in  the  form  of 


tfe-EE 

^  m  n 

(Vym) 


(V;)  (ym-y)  £(vym) 


Ax  Ay 


M’pq  =EE  (*»'  S')(ym'  ?)  £(xn’ym)i=!  4y 
m  n 

(=tl-yL).er 


x*  =  x  ,  y‘  =  y 


=  (x„-  *  ■fx^x„)P(ym"y+y'm  "ymf  f(Vym)Axiy 

m  n  '  •  '  ' 


(x1  ,y‘  )  e  T 
n  m 


i 

1 
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i 


» 


< 


i 


p  q  /  x'  -x  \P  /  y’  -y  5 

S'V»  Vi>  (»^) 


f  (X  y  )ax  Ay 
n  m. 


for  small  variations  of  x*  from  x  and  y1  from  y  we  will  have 

n  m  m 


Mpq  ~  Mpq  +22 


m  n  L 


x'  -  x  y*  -y 
m  n  ,  mm 

p- — = —  +  q 


x  -  x 
n 


y  -y 

m  7 


(x  -x)P(y  -y)  f(x  ,y  )AxAy 
n  m  n  m 


^o'o  =  Moo 


So 


£±a  +  i 
2 

(M'o<J 


_!-E3 -  +  R 

(U00  ) 


where 


^Mqo  ^ 


f  (x  ,  y) Ax-Ay 
n  m 


The  moments  are  also  approximately  invariant  under 
digital  rotation  withe  variation  dependent  upon  the  size  of 
rotation  and  the  interpolation  method  used. 
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4.3.5  Experimental  Results  for  Scene  Analysis 

To  demonstrate  the  advantages  and  limitations  of  moment 
invariants  for  scene  analysis  the  following  experiment  was 
performed.  A  general  segment  of  an  aerial  scene  was  digi¬ 
tized  and  imbedded  into  a  zero  array  since  the  moment  compu¬ 
tation  depends  upon  the  selected  frame  size  which  must  be 
finite.  This  image  was  then  reduced  in  size  by  a  factor  of 
two  by  averaging  four  points  into  one;  inverted  in  magnitude 
modulo  255;  mirror  imaged,  rotated  45°  and  2°  using  nearest 
neighbor  interpolation.  These  images  are  shown  in  Figure  47. 

The  corresponding  invariant  measurements  are  shown  in 
Figure  46.  The  logarithms  have  been  taken  of  the  invariants 
to  compress  the  dynamic  range.  Note  that  the  invariants  are 
close  but  not  identical  for  the  different  transformations  of 
the  original  scene. 

The  main  conclusion  which  can  be  deduced  from  these  ex¬ 
perimental  results  is  the  great  conformity  of  the  theory  with 
the  empirical  results.  This  suggests  moment  invariants  as  a 
practical  and  convenient  tool  in  shape  representation  and 
image  registration.  At  the  end  point  of  the  calculation  one 
obtains  a  set  of  measurements  which  are  invariant  to  any 
translation,  rotation  or  size  change  which  might  occur  in  any 
picture  of  a  scene  due  to  various  factors.  The  method  may 
be  expanded  beyond  this  extent,  i.e.  for  order  p  moments, 
one  has  p+1  invariant  measurements  (refer  to  Section  4.1). 
Hence,  a  large  supply  of  these  invariant  measurements  are 


available . 


EXPERIMENTAL  RESULTS: 

In  this  part  we  verified  empirically  the  seven  invariant 
moments  which  have  been  derived  from  a  given  scene.  The 
measurements,  shown  in  Fig.  4  6  correspond  to  the  pictures  in 
the  following  pages  respectively. 


FINAL  INVARIANTS- 
6.24993 
17.18015 
22.65516 
22.91954 
45.74918 
31 . 83071 
45.58951 


Optical  Window 


FINAL  INVARIANTS- 
6.21612 
16.58215 
21.79714 
22.68471 
44.92849 
31.22309 
44.90914 


Optical  Inverse 


FINAL  INVARIANTS- 
6.31823 
16.80396 
19.72426 
20.43774 
40.52568 
29.31589 
40.47074 


FINAL  INVARIANTS- 
6.22637 
16.95439 
23 . 53142 
24.23687 
48.34990 
32.91619 
48.34356 


Optical  Reduced  2 


FINAL  INVARIANTS- 
6.91980 
19.95532 
26.63924 
26.90140 
53.72453 
37.13457 
53.59021 


Optical  Mirror 


FINAL  INVAR IANTS- 
6.25346 
17.27091 
22.83652 
23.13025 
46.13627 
32.06803 
46.01707 


Optical  Rotation  45° 


Optical  Rotation 


Figure  46 


Optical  Rotation  2  Optical  Rotation 


FIGURE  48. 
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SECTION  V.  HIERARCHICAL  SEARCH  TECHNIQUES 

5.1  Introduction 

The  problem  of  matching  two  images  of  the  same  scene 
taken  by  different  sensors  under  different  conditions  is  a 
challenging  problem  in  scene  analysis.  The  two  scenes  called 
the  reference  scene  and  sensed  scene  may  be  transformed  so 
drastically  by  different  data  collection  geometries  or  sen¬ 
sors  that  it  is  not  possible  to  match  the  scenes.  However, 
in  most  interesting  cases,  the  scenes  can  be  matched  by  a 
photointerpreter.  Only  this  class  of  scenes  will  be  consi¬ 
dered.  A  general  approach  to  matching  two  scenes  involves 
locating  corresponding  regions  in  the  scenes.  If  several 
corresponding  regions  can  be  located,  then  geometric  and  in¬ 
tensity  mappings  may  be  developed  to  match  the  scenes.  There¬ 
fore,  a  basic  problem  is:  given  a  reqion  of  the  reference 
(sensed)  image,  determine  its  location  in  the  sensed  (refe¬ 
rence)  image  (Fig.  48) .  The  standard  computer  approach  to 
this  problem  is  to  use  the  gray  levels  or  a  derived  image  of 
the  reference  region  as  a  template  and  select  the  position  of 
maximum  cross  correlation  between  the  template  at  each  pos¬ 
sible  shift  position  of  the  sensed  image  as  the  match  loca¬ 
tion.  Since  a  template  of  size  M  x  M  can  be  shifted  into 

2 

(N-M)  possible  positions  in  an  N  x  N  image  as  shown  in  Fig¬ 
ure  1,  the  number  of  correlation  computations  can  be  extre¬ 
mely  large.  Fast  correlation  and  edge  correlation  techniques 


N-M 
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FIGURE  43.  TEMPLATE  MATCHING  POSTITIONS  OF  REFERENCE  REGION  OF  SIZE  M  x  M  IN  SENSED  REGION  OF  SIZE  N 
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decrease  the  correlation  time  at  each  shift  position:  however, 

these  methods  still  require  a  computation  at  each  of  the  (N- 
2 

M)  positions.  The  purpose  of  this  section  is  to  describe 
hierarchical  search  techniques  which  are  logarithmetically 
efficient,  i.e.  reduce  the  number  of  search  positions  to  K 
log  (N-M) . 

Although  no  previous  work  on  hierarchical  map  matching 
has  been  discovered,  several  previous  experiments  on  hier¬ 
archical  decompositions  of  pictures  have  motivated  this  study. 

Klinger  and  Oyer[l]  developed  a  regular  decomposition  of 
picture  areas  into  successively  smaller  quadrants  resulting 
in  a  logarithmetic  search.  Tanimoto  and  Pavlidis  [2]  also 
considered  a  hierarchical  data  structure  for  picture  proces¬ 
sing  in  order  to  speed  up  edge  detection  operations. 

Ramaperian  [3]  also  used  multilevel  search  techniques  for 
edge  detection. 

Several  advantages  of  a  structured  approach  are  apparent. 

It  is  not  necessary  to  examine  each  pixel  in  a  high  resolution 
image  to  locate  a  region  at  high  resolution.  The  selecrivity 
of  the  hierarchical  techniques,  especially  coarse-fine  search 
methods,  is  similar  to  the  perceptual  operation  of  the  effi¬ 
cient  human  visual  system.  The  proposed  method  in  which  a 
match  region  is  obtained  at  different  levels  is  extendable  to 
other  problems  such  as  edge  or  object  location.  The  method 
provides  an  arbitrary  degree  of  precision  in  locating  a  region 
which  is  limited  only  by  the  highest  resolution  size  and  the 
uniqueness  of  the  match  region.  Finally,  the  method  permits 

_  ill 
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) 

an  efficient  decomposition  of  the  sensed  scene  into  "infor¬ 
mative"  and  "irrelevant"  regions 


II.  Coarse-Fine  Search  Technique 

In  the  hierarchical  technique  a  structured  set  of 
pictures  at  different  resolutions  will  be  used.  The  high 
resolution  sensed  scene  will  be  denoted,  P  (i,j)  where  for 

1j 

simplicity 

L 

i,j  =  0,1,2,. ..,2  —  1 

The  index  L  is  called  the  level  of  the  search.  The  reference 
region  at  level  L,  QL(i,j)  will  be  assumed  to  be  smaller  than 
the  sensed  scene.  The  agglomerative  rule  by  which  the 
level  K  scene  is  reduced  to  a  level  K-l  scene  is  simple  four 
point  averaging,  i.e. 

Pj^O.j)  =_L_PK'2i,2j)  +  PKf2i,  2j+l)  +PK(2i+l,2i) 

4 

+  PK(2i+l,2j+t) 

K- 

for  i,  j  =  0, 1,  2,  . . .,  2  -  l 

The  reference  image  region  size  varies  with  the  level  and  is 
an  important  consideration.  Obviously,  objects  present  at 
one  level  may  not  be  recognizable  at  a  lower  level.  However, 
allowing  the  reference  image  size  to  change  alleviates  this 
problem.  When  looking  for  a  forest  one  need  not  look  at  the 
leaves.  A  reference  match  region  must  be  selected  for  each 
level  and  the  performance  of  the  algorithm  depends  upon  the 
uniqueness  of  these  reference  regions. 
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A  matching  rule  must  aiso  be  specified  to  guide  the 
search  from  level  K-l  to  level  K.  Several  possible  match 
function  such  as  correlation  may  be  used.  The  match  fun- 
tion  X/^y,  for  x  and  y  should  have  certain  properties, 
such  as 

1.  Identity,  x~x 

2.  Symmetry,  if  x  — y  then  y~x 

3.  Shift  invariance,  if  x-^y  then  x+a  y+a  where  a  is  a 
constant  vector. 

4.  Scale  invariance,  if  x^.y  then  Ay  where  A  is  a  scalar. 

5.  Rotation  invariance,  if  x~-y  then  T  x^y  where  T  is 
a  transformation  matrix. 

6.  Sensor  invariance,  if  x^v  then  t{x)r^y  where  t  is 
the  sensor  transformation. 

The  vectors  x  and  y  may  be  the  gray  level  values  or  derived 
measurements  from  corresponding  sensed  and  reference  regions. 
Note  that  normalized  correlation  satisfies  properties  1-4 
but  not  5  and  6.  Invariant  measurements  such  as  moment  in¬ 
variants  used  with  normalized  correlation  satisfy  properties 
1-6.  Certain  derived  functions  such  as  edges  or  sensor  cor¬ 
rected  images  used  with  invariant  measurements  and  norma¬ 
lized  correlation  may  satisfy  all  the  properties.  Also,  note 
that  other  properties  such  as  invariance  to  weather  or  even 
man  made  changes  may  be  desirable. 

A  specific  example  which  illustrates  a  coarse-fine 
search  will  now  be  considered.  A  quadrant  tree  will  be  used 
for  decomposition  of  the  sensed  scene  as  shown  in  Figure  49. 
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The  efficiency  of  a  quadrant  subdivision  algorithm  for  a 
coarse-fine  search  is  illustrated  in  Figure  50.  Location  of 
an  element  in  an  N  x  N  array,  requires  only  log2 N  steps. 

One  of  the  prominent  features  in  cross-correlation  graphs 
is  the  differences  in  the  location  of  the  peak  when  one  ob¬ 
tains  cross-correlation  of  radar  picture  against  the  corre¬ 
sponding  optical  picture,  and  the  cross-correlation  of  the 
optical  picture  against  the  corresponding  radar  picture.  This 
phenomenon  can  be  explained  fully  by  noting  the  existing 
misregistration  of  the  two  images. 

Consider  two  functions  f(x,y)  and  h(x,y).  The  correlation 
function  of  these  two  can  be  defined  as: 

op 

C(a  ,/?  )  =  ff  f(x  y)  h(x  +  a  .  y  +  0  )  dx  dy 

—  oo 

Now  suppose  there  is  another  function  h(x,y),  which  is  the 
translated  version  of  h(x,y).  Then: 

h ' (x, y)  =  h (x-a, y-b) 

where  a,  b  are  relative  distances  that  has  been  translated. 

The  value  for  correlation  in  this  case  would  be: 


C‘(a  ./?) 
C'(  a  ,/3) 


If1  «*  .  y)  h’{x,  y)  dx  dy 

~jj«x  • y )  h(x  +  a  -  a,  y  +/? -  b)  dx  dy 


— oo 


C(  a  -a  ,  fj  -  b) 


FIGURE  50.  EXAMPLE  OF  QUADRANT  SUBDIVISION  FOR  COARSE-FINE  SEARCH  FOR  DIGITAL  PICTURE. 
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Hence  one  concludes  that  the  new  correlation  function  is  just 
the  translated  version  of  the  original  correlation.  The 
radial  dislocation  is  then 

d  =  /a2  +  b2 

Now  if  f'(x,y)  is  the  translated  version  of  f(x,y),  i.e. 


f ' (x,y)  =  f (x-a  ' ,  y-b ' ! 


then 


OO 

C”  ia.P)  y-b')  h(x+a-a)  (y+y3-b)  dx  dy 


x  =  a '  =  x '  ,  y  -  b  =  y ' 


then 


x  =  x'  +  a'  ,  y  =  y*  +  b' 


C"[OL.fi)  =yy*f(x'.  y')  h{x'+a'-a+a.  y'+b’ -b+0  )  dx'  dy* 


C "(Ct.fi)  =  C  (ot+a'-a,  yS+b'-b) 


Hence  the  relationship  of  this  new  correxation  when  both 
f(x,y)  and  h(x,y)  are  translated  is  just  the  translated 
version  of  the  original  correlation. 


Conclusion: 

In  the  situation  of  the  correlation  graphs  that  was 
referred  before  we  have 

a’  =  -a,  and  b'  =  -b  1 
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hence 


C1  (a  •  (3)  -  C  >9-2b) 


then  the  radial  distance  of  C  with  respect  to  C"  is 


5.2  EXPERIMENTAL  RESULTS 

A  sequence  of  experiments  were  performed  to  investigate 
the  zoom  technique  for  actual  images.  The  simple  4  for  1 
averaging  was  used  for  the  agglomerative  step.  An  inves¬ 
tigation  was  made  of  the  size  reduction  allowable  for  matching 
of  the  large  and  small  stadium  regions  for  both  the  optical 
and  radar  images.  Correlation  was  used  as  the  measure  of 
similarity  and  a  resolution  which  permitted  a  significant 
correlation  peak  was  taken  as  adequate.  Another  experiment 
using  a  derived  edge  picture  and  moment  invariants  was  also 
performed.  Finally,  an  example  is  presented  of  selecting 
moment  measurement  from  the  optical  regions  for  location  of 
the  match  position  in  the  radar  image  at  various  resolutions. 

The  first  sequence  of  experiments  were  performed  to 
verify  the  basic  concept  of  the  zoom  technique.  The  geome¬ 
trically  corrected  optical  image  scanned  at  256  x  256  reso¬ 
lution  was  reduced  to  128  x  128,  64  x  64  and  32  x  32.  A 
similar  sequence  was  produced  for  the  intensity  corrected 
radar  image.  The  following  experiments  were  done: 

Step  1: 

A  64  x  64  optical  picture  was  divided  into  four  equal 
major  regions.  Then  the  invariant  moments  of  the  picture 
were  used  as  a  means  of  representing  the  oicture.  The  top 
right  window  with  (1,1)  coordinates  was  then  used  as  a 
referenced  region  and  the  correlation  function  of  this 
sensed  region  with  the  whole  picture  was  obtained.  The 
reference  region  was  moved  pixel  by  pixel  across  the  whole 
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sensed  picture  and  thus  the  shown  correlation  function  was 
obtained.  The  window  size  was  32  x  32. 

Step  2: 

The  selected  region  in  step  one  was  enlarged  to  64  x  64 
then  this  picture  again  was  divided  into  four  equal  regions. 
The  reigon  with  the  coordinate  (33,  33)  was  selected  and, 
using  a  similar  procedure  as  in  step  one,  the  correlation 
funciton  was  obtained. 

Step  3: 

A  similar  procedure  as  in  step  2  was  performed,  with 
the  exception  of  the  selection  of  the  region,  which  in  this 
case  was  with  the  coordinate  (1,1).  The  obtained  correla¬ 
tion  function  is  shown. 

In  all  the  pictures  shown,  the  photograph  corresponds 
to  the  next  largest  square,  which  contains  the  shaded  region 
For  example,  in  the  Figure  51,  the  photograph  corresponds  to 
the  square  which  is  labeled  c'  in  the  large  square. 

The  resulting  correlations  for  the  optical  image  are 
shown  in  Figures  51,  52  and  53.  The  corresponding  experi¬ 
ment  for  the  radar  images  are  shown  in  Figures  54,  55  and  56 


Selected  Image  Region 
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SECTION  VI.  SCENE  CONTENT  MEASUREMENTS  AND  RECOMMENDA¬ 

TIONS  FOR  FUTURE  WORK 

The  general  problem  of  map  matching  of  images  of  the  same 
scene  taken  under  different  conditions  with  different  sen¬ 
sors  was  considered  in  detail  during  this  study.  Several 
innovative  approaches  and  concepts  were  investigated  and  de¬ 
monstrated  by  experimental  computations  on  images  from  optical 
and  radar  sensors.  Among  these  were: 

*  Techniques  of  geometric  and  sensor  intensity  correc¬ 
tions 

*  Hierarchical  search  for  image  registration 

•  Extraction  of  invariant  features  and  invariant  feature 
correlations 

•  Edge  extraction  and  edge  correlations 

Since  most  correlation  techniques  work  well  on  locating 
a  region  from  an  image,  one  approach  to  map  matching  of  two 
different  images  is  to  develop  inverse  transformations  to  make 
the  images  appear  as  a  region  from  one  image.  This  general 
transformation  depends  upon  both  sensor  and  imaging  characte¬ 
ristics  and  is  difficult  if  not  impossible  to  determine..  How¬ 
ever,  it  was  possible  to  develop  simple  approximations  which 
rectify  not  only  the  geometric  but  also  the  intensity  of 
two  arbitrary  images.  A  "map  warping"  approach  which  involves 
a  two  dimensional  polynomial  mapping  from  corresponding 
points  in  the  two  images  was  used  for  geometric  corrections. 

A  new  technique  based  upon  the  Karhunen-Loeve  transformation 
was  used  to  determine  a  principal  component  of  the  image  with 
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a  corrected  intensity  distribution.  This  approach  is  an 
example  of  a  class  of  intensity  matching  transformations 
which  because  of  their  effectiveness  and  simplicity  may  be¬ 
come  important  techniques  for  map  matching.  These  techni¬ 
ques  only  require  that  intensity  statistics  be  assumed  or 
experimentally  determined.  Other  information  that  could  be 
obtained  about  the  scene,  such  as  object  composition  or 
height,  would  permit  the  develop  of  more  sophisticated  methods. 

Hierarchical  search  techniques  for  image  registration 
were  introduced  in  this  report  and  promise  to  be  a  highly 
efficient  and  effective  method  for  map  matching.  An  impres¬ 
sive  feature  of  these  methods  is  logarithmic  computational 
efficiency.  A  special  case  of  the  hierarchical  techniques 
permits  a  "zoom"  to  be  used  to  locate  an  object  of  interest. 
That  is,  an  image  can  be  searched  at  low  resolutions  and  a 
region  is  located  approximately.  Then  the  approximation  can 
be  refined  to  any  desired  level  by  investigating  selected 
regions  in  higher  resolution  images.  The  hierarchical  tech¬ 
nique  was  demonstrated  for  locating  regions  from  an  optical 
image  in  an  optical  image,  regions  in  a  radar  image  in  a 
radar  and  regions  of  a  radar  image  in  an  optical  image.  The 
matching  of  radar  and  optical  image  requires  the  extraction 
of  measurements  which  are  invariant  to  sensor  and  geometric 
distortions.  This  requirement  led  to  the  study  of  invariant 
measurements. 
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Certain  measurements  have  a  high  degree  of  invariance 
to  the  types  of  geometric  and  sensor  transformations  produced 
in  optical  and  side  looking  radar  imaging.  Moment  measure¬ 
ments  of  the  image  intensity  can  be  transformed  into  a  set  of 
invariants  which  are  invariant  to  translation,  rotation,  scale 
and  other  transformations.  These  invariants  were  computed 
and  used  as  measurements  in  correlation  experiments.  It  was 
noted  that  due  to  the  digital  nature  of  the  computation  the 
measurements  were  only  approximately  invariant  but  did  work 
quite  well  for  map  matching  optical  and  radar  images. 

Edge  structure  also  remains  invariant  to  many  sensor 
transformations  and  thus  edge  correlation  for  map  matching 
was  also  considered  in  detail.  A  new  decomposition  method 
for  3  by  3  windows  was  developed  and  used  for  edge  detection. 
Although  edge  structure  is  very  invariant  to  sensor  trans¬ 
formations,  the  edge  correlator  may  be  sensitive  to  geometric 
transformations  depending  uoon  the  edge  detection  technique. 
Experimental  optical-radar  correlations  were  performed  to 
demonstrate  the  efficiency  of  the  edge  correlation  approach. 

A  common  problem  encountered  in  all  correlation  methods 
is  a  sensitivity  of  performance  to  the  content  of  the  scene. 
This  phenomena  was  demonstrated  in  several  examples  and 
requires  further  investigation. 

Each  of  these  techniques  and  developments  have  been 


successfully  applied  to  a  limited  set  of  images.  The  use  of 
a  limited  number  of  images  was  necessary  at  this  time  in  order 
to  provide  a  well-defined  environment  for  the  development 
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of  the  basic  techniques. 

To  demonstrate  the  true  usefulness  of  these  approaches, 

they  must  be  applied  to  many  scenes  of  interest.  Therefore, 

\ 

it  is  recommended  that  the  following  work  be  performed: 

1.  Establish  the  performance  of  these  approaches  to 
sequential  map  matching  using  hierarchical  search, 
edge  features,  and  structured  invariant  moments  with 
respect  to  terrain  and  type. 

2.  Identify  the  types  of  scenes  and  features  which  are 
best  suited  for  these  approaches. 

3.  Perform  theoretical  and  experimental  statistical 
analysis  on  these  approaches  and  determine  suita¬ 
bility  of  these  approaches  to  the  perturbed  scenes. 

4.  Verify  the  validity  of  experimental  performances 
through  statistical  analysis . 


